{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Boundary conditions tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial aims to demonstrate how users can implement various boundary conditions in Devito, building on concepts introduced in previous tutorials. Over the course of this notebook we will go over the implementation of both free surface boundary conditions and perfectly-matched layers (PMLs) in the context of the first-order acoustic wave equation. This tutorial is based on a simplified version of the method outlined in Liu and Tao's 1997 paper (https://doi.org/10.1121/1.419657).\n",
    "\n",
    "We will set up our domain with PMLs along the left, right, and bottom edges, and free surface boundaries at the top as shown below.\n",
    "\n",
    "<img src=\"figures/boundary_conditions.png\" style=\"width: 220px;\"/>\n",
    "\n",
    "Note that whilst in practice we would want the PML tapers to overlap in the corners, this requires additional subdomains. As such, they are omitted for simplicity.\n",
    "\n",
    "As always, we will begin by specifying some parameters for our `Grid`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "shape = (101, 101)\n",
    "extent = (2000., 2000.)\n",
    "nbpml = 10  # Number of PMLs on each side"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will need to use subdomains to accomodate the modified equations in the PML regions. As `Grid` objects cannot have subdomains added retroactively, we must define our subdomains beforehand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from devito import SubDomain\n",
    "\n",
    "class MainDomain(SubDomain):  # Main section with no damping\n",
    "    name = 'main'\n",
    "    def __init__(self, PMLS):\n",
    "        super().__init__()\n",
    "        self.PMLS = PMLS\n",
    "            \n",
    "\n",
    "    def define(self, dimensions):\n",
    "        x, y = dimensions\n",
    "        return {x: ('middle', self.PMLS, self.PMLS), y: ('middle', 0, self.PMLS)}\n",
    "\n",
    "\n",
    "class Left(SubDomain):  # Left PML region\n",
    "    name = 'left'\n",
    "    def __init__(self, PMLS):\n",
    "        super().__init__()\n",
    "        self.PMLS = PMLS\n",
    "\n",
    "    def define(self, dimensions):\n",
    "        x, y = dimensions\n",
    "        return {x: ('left', self.PMLS), y: y}\n",
    "\n",
    "\n",
    "class Right(SubDomain):  # Right PML region\n",
    "    name = 'right'\n",
    "    def __init__(self, PMLS):\n",
    "        super().__init__()\n",
    "        self.PMLS = PMLS\n",
    "\n",
    "    def define(self, dimensions):\n",
    "        x, y = dimensions\n",
    "        return {x: ('right', self.PMLS), y: y}\n",
    "    \n",
    "    \n",
    "class Base(SubDomain):  # Base PML region\n",
    "    name = 'base'\n",
    "    def __init__(self, PMLS):\n",
    "        super().__init__()\n",
    "        self.PMLS = PMLS\n",
    "\n",
    "    def define(self, dimensions):\n",
    "        x, y = dimensions\n",
    "        return {x: ('middle', self.PMLS, self.PMLS), y: ('right', self.PMLS)}\n",
    "\n",
    "\n",
    "main_domain = MainDomain(nbpml)\n",
    "left = Left(nbpml)\n",
    "right = Right(nbpml)\n",
    "base = Base(nbpml)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And create the grid, attaching our subdomains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from devito import Grid\n",
    "\n",
    "grid = Grid(shape=shape, extent=extent,\n",
    "            subdomains=(main_domain, left, right, base))\n",
    "x, y = grid.dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then begin to specify our problem starting with some parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "density = 1.  # 1000kg/m^3\n",
    "velocity = 4.  # km/s\n",
    "gamma = 0.0002  # Absorption coefficient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need a `TimeFunction` object for each of our wavefields. As particle velocity is a vector, we will choose a `VectorTimeFunction` object to encapsulate it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from devito import TimeFunction, VectorTimeFunction, NODE\n",
    "\n",
    "p = TimeFunction(name='p', grid=grid, time_order=1,\n",
    "                 space_order=6, staggered=NODE)\n",
    "v = VectorTimeFunction(name='v', grid=grid, time_order=1,\n",
    "                       space_order=6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `VectorTimeFunction` is near identical in function to a standard `TimeFunction`, albeit with a field for each grid dimension. The fields associated with each component can be accessed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  ...\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]]\n",
      "\n",
      " [[0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  ...\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]\n",
      "  [0. 0. 0. ... 0. 0. 0.]]]\n"
     ]
    }
   ],
   "source": [
    "print(v[0].data)  # Print the data attribute associated with the x component of v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may have also noticed the keyword `staggered` in the arguments when we created these functions. As one might expect, these are used for specifying where derivatives should be evaluated relative to the grid, as required for implementing formulations such as the first-order acoustic wave equation or P-SV elastic. Passing a function `staggered=NODE` specifies that its derivatives should be evaluated at the node. One can also pass `staggered=x` or `staggered=y` to stagger the grid by half a spacing in those respective directions. Additionally, a tuple of dimensions can be passed to stagger in multiple directions (e.g. `staggered=(x, y)`). `VectorTimeFunction` objects have their associated grids staggered by default.\n",
    "\n",
    "We will also need to define a field for integrating pressure over time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_i = TimeFunction(name='p_i', grid=grid, time_order=1,\n",
    "                   space_order=1, staggered=NODE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we prepare the source term:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from examples.seismic import TimeAxis, RickerSource\n",
    "\n",
    "t0 = 0.  # Simulation starts at t=0\n",
    "tn = 400.  # Simulation length in ms\n",
    "dt = 1e2*(1. / np.sqrt(2.)) / 60.  # Time step\n",
    "\n",
    "time_range = TimeAxis(start=t0, stop=tn, step=dt)\n",
    "\n",
    "f0 = 0.02\n",
    "src = RickerSource(name='src', grid=grid, f0=f0,\n",
    "                   npoint=1, time_range=time_range)\n",
    "\n",
    "# Position source centrally in all dimensions\n",
    "src.coordinates.data[0, :] = 1000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "tags": [
     "nbval-skip"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAF9CAYAAACH0lvIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nO3deZxkZXn3/89VvW+zLzgwW7PJgApmRBQVRCNGfSB5SDQqCiaKJsZoEhLlZx5jEA1u6M9oDJjEDWOMBB+IGkT2iKAZ2WTYZoZZGAaYnq3X6q7q6uv545zTU1NTXX26u9Yz3/frVa/qPnXOqftM9VRddd3Xfd/m7oiIiIhMJ1XrBoiIiEhjUNAgIiIisShoEBERkVgUNIiIiEgsChpEREQkFgUNIiIiEktzrRtQ75YsWeJr1qypdTNERESq4le/+tUed19a7DEFDdNYs2YNGzZsqHUzREREqsLMtk/1mLonREREJBYFDSIiIhKLggYRERGJRUGDiIiIxKKgQURERGJR0CAiIiKxNETQYGbHmNnfm9k9ZjZiZm5ma2IemzKzy8xsm5mNmtmDZnZBZVssIiKSPA0RNADHAW8G9gP/PcNjPwF8HPgy8FvAvcD3zewN5WygiIhI0jXK5E53uftyADN7N/C6OAeZ2TLgUuBKd/9cuPl2MzsOuBL4cSUaKyIikkQNkWlw94lZHnou0ApcW7D9WuAFZrZ2Tg0TERE5gjRE0DAHJwNjwOaC7RvD+3XVbY6IiEjjSnrQsAg44O5esH1f3uMyQ7kJ5+Gn+8lNFP6ziohIkiU9aJgVM7vEzDaY2Ya+vr5aN6fu/OD+p3nT3/+Mcz5/B7sOpGvdHBERqZKkBw37gQVmZgXbowzDPopw92vcfb27r1+6tOjqoEe0+3bsp6XJ2L53hFsffa7WzRERkSpJetCwEWgDji3YHtUyPFLd5iTDxl0DrF+9iJ62Zp54bqjWzRERkSpJetBwE5AF3l6w/ULgYXffWv0mNbbx3ASPPTPAySvmcfzybjbtHqx1k0REpEoaZZ4GzOx3wx9/I7z/LTPrA/rc/c5wn3Hgm+7+hwDuvtvMrgIuM7NB4D7gLcA5wHlVvYCE2NI3zNj4BCcfPY/B0XFuUfeEiMgRo2GCBuD7Bb//Q3h/J3B2+HNTeMv3UWAI+CBwFPA48GZ3/2FlmplsG3f1A3DKivnsHcrwvQ1PsXdojMXdbTVumYiIVFrDBA3uXljMGGsfd88BV4Q3maNHnxmgrTlF79JunukfBWDT7iEFDSIiR4Ck1zRIme3qH2XFgg6aUsbxy7sB2PSc6hpERI4EChpkRvoGx1gaZhWOmtdOW3OKHftGatwqERGpBgUNMiN7BsdYOi8IGsyMJd1t7B3K1LhVIiJSDQoaZEbyMw0AS3ra6Bsaq2GLRESkWhQ0SGzpTI7BsXGWzTsYNCztbmWPMg0iIkcEBQ0SW99gkFHIzzQs7mpjrzINIiJHBAUNElvfUDDEcmlPfvdEK3uHM0xoxUsRkcRT0CCx7R4IMgrLetonty3pbiM34RxIZ2vVLBERqRIFDRJbVPCYn2mIJnVSF4WISPIpaJDY+gbHaEoZi7paJ7ct6Q5+1ggKEZHkU9AgsfUNjrG4q5Wm1MHZuqOiSI2gEBFJPgUNEtvuwbFDuiZA3RMiIkcSBQ0SW7HVLBd0tNCUMvYoaBARSTwFDRLbwOg4CzpaDtmWShmLu1rZM6juCRGRpFPQILH1p7PM6zh8NfXF3W3sHVamQUQk6RQ0SCzuTn86y/yCTAMEXRQD6fEatEpERKpJQYPEMpzJkZvwokHDvI5m+jW5k4hI4ilokFgGwqBgXnuRoKG9hYFRBQ0iIkmnoEFiiTIJxTMNLZNBhYiIJJeCBomlZNDQ3sJwJsd4bqLazRIRkSpS0CCxREHDvClqGgAGR1UMKSKSZAoaJJaBaTINgOoaREQSTkGDxFI60xAGDRp2KSKSaAoaJJaBdBYz6Gk7fHKnee3BNmUaRESSTUGDxNKfzjKvvYVU3gqXkfmdUaZBQYOISJIpaJBYBkbHi04hDappEBE5UihokFimmkIaVNMgInKkUNAgsZQKGrpam0gZmkpaRCThFDRILKWCBjMLZoVU94SISKI1RNBgZivN7Doz6zezATO73sxWxTx2lZl908x2mFnazJ4wsyvMrKvS7U6SgbAQcirz2jWVtIhI0hWvbKsjZtYJ3AaMARcBDlwB3G5mL3T34RLHdgG3AC3A/wF2AC8B/hY4HnhLZVufHAOj2aJzNETmdTQzoBkhRUQSre6DBuA9QC9wortvBjCzh4BNwHuBq0oceyZBcHCuu98cbrvdzBYBl5pZp7uPVK7pyZDNTTCanaC7yBwNEWUaRESSrxG6J84D7o0CBgB33wrcDZw/zbGt4f1AwfYDBNd++KQDcpjhsSCDUCpomK+aBhGRxGuEoOFk4OEi2zcC66Y59haCjMSnzWydmXWb2TnAB4F/LNW1IQcNxQgaetqbNXpCRCThGiFoWATsL7J9H7Cw1IHuPgq8guA6NwKDwK3AD4E/meo4M7vEzDaY2Ya+vr7ZtjsxhsdyAHSVCBq621om9xMRkWRqhKBh1sysHfgesAx4B3AW8JcEBZBfmeo4d7/G3de7+/qlS5dWpa31bGgsyCB0t08dNHS1NTGcGcfdq9UsERGpskYohNxP8YzCVBmIfH8InA0c5+5bwm13mVk/cI2Z/aO7P1i2libUUJhB6G5rmnKfrrZm3GEkkyuZkRARkcbVCJmGjQR1DYXWAY9Mc+wLgP15AUPkl+H9SXNs2xEhKoQsFQxEj0X7iohI8jRC0HAjcIaZ9UYbzGwNwXDKG6c59llgoZkdV7D9peH902VqY6INhfMvdLWWqmkIshBDChpERBKrEYKGrwHbgBvM7HwzOw+4AXgKuDraycxWm9m4mX0s79hvEBQ//tjMLjKzV5vZXwKfA35FMGxTphEFAj2lahpao0yDiiFFRJKq7oOGcFjkOcATwLeB7wBbgXPcfShvVwOayLsmd98GnAE8QDCL5I8JJou6BvhNd5+owiU0vDjdE9FwTGUaRESSqyEq1tx9B3DBNPtso8hkTe7+CPDmyrTsyDCUGae1OUVL09QxZhRQjGQUNIiIJFXdZxqk9oZGx+mZZkRElzINIiKJp6BBpjU8Nj7tMMruNtU0iIgknYIGmdbQ2PRzL3SFoyc05FJEJLkUNMi0hsayJSd2goOjJ9Q9ISKSXAoaZFrDY7mSi1UBpFJGZ2uTMg0iIgmmoEGmFaemAYJiyGGNnhARSSwFDTKtwbHxaTMNEBRDDqkQUkQksRQ0yLSGYwYNXW3qnhARSTIFDVJSbsJjr1zZ1dqsQkgRkQRT0CAlRTUKcbsnlGkQEUkuBQ1SUpx1JyJdChpERBJNQYOUdDBoKD1PQ7CPCiFFRJJMQYOUFE0LHU3eVEqX5mkQEUk0BQ1S0kgmCBo6Y2Ya0tkcuQmvdLNERKQGFDRISdFS150xMg2Ti1ZpgicRkURS0CAlRZmGrtZ4mQaAEdU1iIgkkoIGKSnKNHTECBo6w31GlGkQEUkkBQ1S0sFMw/TdEx2TQYMyDSIiSaSgQUqaSSFkR0uwTzqroEFEJIkUNEhJI5lxmlJGa9P0fyqdyjSIiCSaggYpaXgsR2drE2Y27b5R90RaNQ0iIomkoEFKSmdykxmE6UTDMpVpEBFJJgUNUtJwZjxWESSoe0JEJOkUNEhJI5lcrOGWkN89oaBBRCSJFDRISSMzyTS0KNMgIpJkChqkpJlkGpqbUrQ2pRjJqhBSRCSJFDRISSOZXKxlsSMdrU3qnhARSSgFDVLSyNg4HS3xuicgKIZU94SISDI1TNBgZivN7Doz6zezATO73sxWzeD4k8zs+2a2x8zSZva4mX2wkm1OgpGsMg0iIhKI/xWyhsysE7gNGAMuAhy4ArjdzF7o7sPTHL8+PP4O4N1AP3A80F3BZifCyFgu1rLYkc7WJk0jLSKSUA0RNADvAXqBE919M4CZPQRsAt4LXDXVgWaWAr4F3Oruv5P30O2Va24yZHMTZHITsSd3AuhsadYqlyIiCdUo3RPnAfdGAQOAu28F7gbOn+bYs4GTKBFYSHGTi1XNIGhQ94SISHI1StBwMvBwke0bgXXTHPuK8L7dzO41s6yZ7TazL5lZR1lbmTDpyaBBhZAiItI4QcMiYH+R7fuAhdMcuyK8/x5wM/CbwGcIahv+tdgBZnaJmW0wsw19fX2za3ECDIfdDDMthFTQICKSTI1S0zAXUWB0rbt/LPz5DjNrAq40s5Pc/dH8A9z9GuAagPXr13v1mlpfRsaCD/+OlhkEDS0qhBQRSapGyTTsp3hGYaoMRL694f1PC7bfHN6fNod2JdrIZKZhpt0TKoQUEUmiRgkaNhLUNRRaBzwS49hSJmbVoiNA1M0QdxrpYN9mRrMTTEwcsQkaEZHEapSg4UbgDDPrjTaY2RrgzPCxUv6LYH6Hcwu2vz6831CeJibPbEZPRPuqi0JEJHkaJWj4GrANuMHMzjez84AbgKeAq6OdzGy1mY2bWVS7gLvvBf4OeJ+ZfcrMXmtmHwE+BnwzfxinHGqye2KGoyeCYxU0iIgkTUMUQrr7sJmdA3wB+DZgwK3Ah9x9KG9XA5o4PBi6HBgE/hi4FHgG+CzwiQo3vaFF2YL2GRZCApqrQUQkgRoiaABw9x3ABdPss40gcCjc7gSTO2mCpxlIz6p7IviT0vLYIiLJ0yjdE1IDs8k0qHtCRCS5FDTIlNKZHG3NKZpShyVvphSNtFD3hIhI8ihokCmls7kZDbcEZRpERJJMQYNMaSSTo3MGXROgIZciIkmmoEGmlM7maJ9hpqEjLIRMa1ZIEZHEUdAgU0pncjMaOQFMZibUPSEikjwKGmRK6UxuRotVwcFCSAUNIiLJo6BBpjSSzU12N8TV1pwiZRo9ISKSRAoaZEqjmRwdLTP7EzEzOlqalGkQEUkgBQ0ypZHs+OQMjzPR0dpMWjNCiogkjoIGmVI6MzGj2SAjna3KNIiIJJGCBplSOjM+40JIUNAgIpJUChqkKHcnnZ35kEsIRlCoEFJEJHkUNEhRY+MTTDgznkYaokyDahpERJJGQYMUNRpOAz2b7omOlmZ1T4iIJJCCBikqWjtitpkGrT0hIpI8ChqkqChTMJuaBhVCiogkk4IGKSoqZJzNkEsVQoqIJJOCBikq6l6YbaYhnc3h7uVuloiI1JCCBikqyhTMbp6GZnITTiY3Ue5miYhIDSlokKKimoTZFEJGgYa6KEREkiV20GCB88zsc2b2dTNbHW4/y8xWVK6JUgtzGXLZqeWxRUQSKdZqRGa2EPgx8FJgEOgG/h7YDrwH2Af8aYXaKDVwcPTEbBasUtAgIpJEcTMNnwVWAmcCiwHLe+wW4DVlbpfUWHpOmYYg0FD3hIhIssT9Gnk+cKm732NmhZ8iOwgCCkmQdDgN9FxqGjSVtIhIssTNNHQDT0/xWDuHZh4kAdLZHE0po6Vp5i/tZPeEZoUUEUmUuEHD48DrpnjsLODX5WmO1IuRTI6OlibMZh40RIWQ6p4QEUmWuN0T/wB82cz6gX8Nty0ws3cBfwJcUonGSe2MZnOz6poAjZ4QEUmqWJkGd78GuAr4W2BzuPmnwDXAF939O5Vp3kFmttLMrjOzfjMbMLPrzWzVLM7zETNzM/tZJdqZFFGmYTY6JjMNqmkQEUmS2OPp3P0jZvZV4DeBZcBe4Kfu/mSlGhcxs07gNmAMuAhw4ArgdjN7obsPxzxPL/DXwO5KtTUp0pncrKaQhoOjJ5RpEBFJlhkNwnf37cA/VagtpbwH6AVOdPfNAGb2ELAJeC9BFiSOrwLfAU5khtd+pElnc7NarAqgM5oRUoWQIiKJMuUH50xT/+6+Y+7NmdJ5wL1RwBA+31Yzu5tgOOi0QYOZvQ14MfBW4PpKNTQp5pJpSKWMtuaUCiFFRBKm1LftbQTdAHHN7hMmnpOBG4ps3wj83nQHhzNafgH4K3ffN5sRAUeadDbH/I6WWR/fEa50KSIiyVEqaPgDDgYNbQS1AAPAvwPPAUcBbwZ6gE9UsI0Ai4D9RbbvAxbGOP6zwBPAN+I8mZldQjgiZNWqGddaJkI6M/vRExBM8KRMg4hIskwZNLj7N6KfzeyLwH3A77i7522/HPi/wLoKtnFOzOyVwDuBF+e3vZRwtMg1AOvXr59JtiUx0tnZj56AINOgyZ1ERJIl7uRObwWuLvzQDX//R+Bt5W5Ygf0UzyhMlYHIdzXwz8BOM1tgZgsIgqWm8Pe28jY1GUbmUNMAQaZhVJkGEZFEiTuCoBtYOsVjy4Cu8jRnShsJ6hoKrQMemebYk8Lb+4o8th/4M+CLc2pdAqWzOdrnEDR0tjZpyKWISMLEDRruAD5lZo+6+/9EG83sdOCT4eOVdCPwOTPrjeaFMLM1BKtufmSaY19dZNsXCQo3P8DByaoklJtwMuMTdLbMflRqe0sTg6Oa3ElEJEnifir8CcES2Pea2VMEhZDLCVa33Bo+XklfC5/jBjP7a4ICzU8ATxF0PwBgZquBLcDl7n45gLvfUXgyMzsANBd7TPKWxW6N23t1uI6WJvoGx8rVJBERqQNxp5HeCjyfIMV/K8FskLcSTKx0krtvq1QDw+cfBs4hGAHxbYIJmrYC57j7UN6uRpBBmP2nnUwuad3ROvtMg7onRESSZybTSGcJvvF/rXLNKfn8O4ALptlnGzGW6Xb3s8vTqmQazUwAzHn0hOZpEBFJFn0jl8OMZINMw9xGTzRrngYRkYSJlWkws62Unh3S3f3Y8jRJai36sJ9bpiFFOpvD3dEMnCIiyRC3e+JODg8aFgMvB4YIVqCUhIiChtkuWAVBwJGbcDK5CdqaKznDuIiIVEusoMHdLy62PZwo6SaCkRWSEFEtwpy6J8IiytGMggYRkaSYU02Dux8gWNfhY+VpjtSDaNTDXNeeAC2PLSKSJOUohBwFjinDeaROTM7TMIfuiShLEQ3fFBGRxjfrgfhm1gycAnycYJpnSYjR7NwzDe3KNIiIJE7c0RMTTD16YgB4Y9laJDUXdU/MpaYhOlbDLkVEkiNupuFyDg8aRoHtwH+5e39ZWyU1NTl6Yg4FjFGWQpkGEZHkiDt64uMVbofUkXQ2R3tLilRq9vMrTBZCKtMgIpIYsQohzew2M3v+FI+dYGaapyFB0pncnIogQZkGEZEkijt64mxg3hSP9QBnlaU1UhdGMjk657BYFSjTICKSRDMZcjlVIeSxBLNCSkKMht0Tc3FwyKWCBhGRpJjy66SZvQt4V/irA9eY2WDBbh0Ewy5vrUzzpBZGMuNzzjRoyKWISPKU+jo5AeTCmxX8Ht32Al8F/rCyzZRqSmfnXtPQ1pzCTN0TIiJJMuXXSXf/JvBNADO7Hfgjd3+sWg2T2klncizobJ3TOcyMzpYmZRpERBIk7pDLV1e6IVI/RjI5ViyY+yJTHa0KGkREkqRUTcM7gR+5+97w55Lc/VtlbZnUzEgZhlxCGDSoe0JEJDFKZRq+AZxBULfwjWnO44CChoQYzebmtO5EpKNFQYOISJKUChrWAs/k/SxHiGCehnJkGpoZUfeEiEhilCqE3F7sZ0m2iQkPRk/MccglQEdLilFlGkREEmNuM/hI4oyOh8til6OmQaMnREQSpVQh5FamngWykLv7seVpktRSOZbFjnS2NjOSGZnzeUREpD6UykHfSfygQRIiKlwsRyFke0sTo9mJOZ9HRETqQ6mahour2A6pE1F3QnkyDU2MZMbnfB4REakPqmmQQ0TdE2Wbp0E1DSIiiRE7aDCz483sm2b2hJkNh/ffMLPjKtlAqa4oM1DO7omJCfVyiYgkQaygwczOBh4E3gTcC/xDeP+/gF+b2VmVaqBU1+hk98Tch1xGXRzRiAwREWlscTMNnwfuB1a7+zvd/S/d/Z3AGuCB8PGKMrOVZnadmfWb2YCZXW9mq2Ict97MrjGzx8xsxMx2mNl3zEwTVhVRztETUReHZoUUEUmGuEHDOuDT7j6Uv9HdB4FPAyeXu2H5zKwTuA14PnAR8A7geOB2M+ua5vDfD9v3JeC3gI8ALwY2mNnKijW6QZW7piH/nCIi0tji5qB3AlOtldwKPF2e5kzpPUAvcKK7bwYws4eATcB7gatKHPtpd+/L32BmdwNbw/N+rCItblDlHHIZBR6jKoYUEUmEuJmGTwN/a2Yr8jea2dHA3wCfKnfDCpwH3BsFDADuvhW4Gzi/1IGFAUO4bTvQBxxd5nY2vHIPuQRlGkREkiJupuEsYB7wpJndCzwHLCdYBfM54OywWBKC2SEvKnM7TwZuKLJ9I/B7Mz2ZmZ0ELAMenWO7Eif6gG9vLmNNgzINIiKJEDdoeAUwTrDq5erwBgdXwXxl3r6VGF+3CNhfZPs+YOFMTmRmzcA/EmQa/nmKfS4BLgFYtWraWstESWfGaW9JkUrZnM/V3qqgQUQkSWIFDe6epJEGXwZeDrzR3YsFIrj7NcA1AOvXrz+iJhkIlsWe+3BLONg9odETIiLJUJ5Ph8rbT/GMwlQZiKLM7EqCDMJF7n5zmdqWKOlsriwjJ0BDLkVEkmZGQUM4RHEl0F74mLvfVq5GFbGR4sM61wGPxDmBmX0U+DDwAXf/dhnblijpTK4sRZCQN+RS3RMiIokQK2gws17gO8Dp0abw3sOfHSjPJ01xNwKfM7Ned38ybNMa4EyCeRdKMrM/Ba4APuruX65gOxveSDmDhmjIpTINIiKJEDfT8E/AKuBDwGNApmItKu5rwJ8AN5jZXxMEKZ8AngKujnYys9XAFuByd7883Pb7wBeBm4DbzOyMvPMOuHusTMWRIp3J0V7m7gkNuRQRSYa4QcNLgIvd/T8q2ZipuPuwmZ0DfAH4NkF241bgQwWzVBpBxiN//onXh9tfH97y3QmcXaFmN6R0NseS7qnm8ZqZ5qYUrU0pjZ4QEUmImcwIWe3swiHcfQdwwTT7bONg10m07WLg4kq1K2lGMuN0tnaW7XztLSnNCCkikhBxZ4T8FPDhGOs8SINLZ3JlmUI60tnaPLnctoiINLa48zR828yeD2wLZ4QsHOZYiVkgpQZGyjjkEoIRFOnsRNnOJyIitRN39MTFwGVAjmCFyMKuiiNqAqQkK+foCQiKIdPKNIiIJELcmoa/BX4A/KG7H6hge6SGchNOZnyirN0TQaZBNQ0iIkkQt6ZhMfAPChiSrZwrXEaCTIOCBhGRJIgbNPwMOKmSDZHaiwoWy13ToHkaRESSIW73xAeBfzez/QSTJB223oO7q9qtwUUZgY4yLVgFYaZB3RMiIokQ99Ph0fD+WyX2qeQ00lIFleie6FSmQUQkMeIGDZejERKJNzKZaShzIaSCBhGRRIg7T8PHp3rMzM4G3lmm9kgNTXZPlLGmoSuc3MndMbPpDxARkboVtxDyEGZ2nJldbmZbCdaAeHN5myW1EGUayjp6orWJCYexcZW8iIg0uthBg5nNN7NLzOxu4HHgowQFkX8MrKhQ+6SKKlXTAKiLQkQkAUoGDWaWMrM3mNn3gGeAfwRWA18Jd/mQu1/t7gMVbqdUQTRzYzlHT0RBw4hGUIiINLwpPx3M7PPA24BlwCjBjJDfBG4B5gF/Uo0GSvWMVKCmIQpARsY0lbSISKMr9ZXyzwhGTPwYuNjd90YPmJlGUiRQJWoauqJMg7onREQaXqnuiX8GBoE3Ao+b2ZfN7PTqNEtqYTSbwwzammdVH1tUh4IGEZHEmPLTwd3fAxwFvB3YALwXuMfMHgU+jOZtSJyRTI7OlqayDo3sDLsn0ll1T4iINLqSXyndfdTdv+vurwdWcXB57I8ABlxpZheaWXvlmyqVNpLJlbUIEvIKIZVpEBFpeLHz0O7+jLt/xt1PAU4nGEFxPMHU0s9UqH1SRenMOB2t5euagLygYUxBg4hIo5vVJ4S7b3D3DxDMz3ABcEc5GyW1EXRPlDvTEI6eyKh7QkSk0c3pE8LdswRDMX9QnuZILaWzubKuOwGap0FEJEnKm4uWhpbO5Mo63BKCkRhmmhFSRCQJFDTIpJFMrqwTOwGYGZ0tWh5bRCQJFDTIpEp0TwB0tjWrpkFEJAEUNMikSnRPQFDXoEyDiEjjU9Agk0Yy45OjHcqpQ90TIiKJoKBBJqWzOdrLXNMAQaZBhZAiIo1PQYMAkM1NkM15RbonulTTICKSCA0TNJjZSjO7zsz6zWzAzK43s1Uxj203s8+a2TNmljaze8zsVZVucyNJZ8u/wmVE3RMiIsnQEEGDmXUCtwHPBy4C3kEwhfXtZtYV4xT/DLwH+BjwJoJpr39iZqdWpsWNJ+o+qMjoCRVCiogkQvmr3irjPUAvcKK7bwYws4eATQSrb1411YFm9iLgbcAfuPvXw213AhuBy4HzKtv0xhB9qFck09DarKBBRCQBGiLTQPDBfm8UMAC4+1bgbuD8GMdmge/lHTsO/Btwrpm1lb+5jWd4LKg56KrA6ImgEFI1DSIija5RMg0nAzcU2b4R+L0Yx25195Eix7YCx4U/V9y2PcP0DY2xvKedZfPaKjJSYbYmg4a28v9JdLU2MZLN4e6YWdnPP1sjmXF27k8zkM4yODrOwGiWsewEOXcm3JmYcHITzoQT/O5e6yaLiBxmcVcbF/zGMVV5rkYJGhYB+4ts3wcsnMOx0eOHMLNLgEsAVq2KVWsZy3W/2smXb59MltC7tItTj1nAGb2LOfeUo5jf0VK255qpSndPuMNodqIiNRNx7dw/wu2P7ebOJ/p4aGc/uwfHatYWEZFyWfe8eQoaasndrwGuAVi/fn3Zvl6+7aWrOH3tIp4dGGXXgTQPPz3AXZv2cP39T/OxGx/mraev4kOvOYH5ndUPHoYzFcw0tIUrXWbGaxI0PNk3xKdveoybH3kOd1i5qINXnbCUtUu6WLmok4WdLXS3NdPT3kx7SxNNKaPJDDM7+HMKUmbUT55ERCSQqmIGt1GChv0UzyhMlUUoPHb1FMfCwYxDxa1Y0MGKBR2HbHN3HtrZz7fu2c637tnODx96hi/9/mm87NjF1WoWACNjQaahEkFDNMvk8FiOxd1lP40gc0MAACAASURBVP2U3J1v/nwbn/rxY7Q1p3j/2cdxwW8cw5rFnXXVTSIi0igapRByI0FtQqF1wCMxjl0bDtssPDYDbD78kOoxM160cgGff/OLuOH9ZzK/o4V3/ssv+OFDu6rajqHJQsjyZwK6w0zDcBWLId2dT/7oUT7+n4/wyuOXcOulZ3HpuSeydkmXAgYRkVlqlKDhRuAMM+uNNpjZGuDM8LFS/hNoIa9g0syagbcAN7t73XRsn3L0fP7jj17OaSsX8mffe4Cfb95TteeOZmysxNoTBzMN1QsavnL7Zv7pZ1u5+OVr+No717Osp71qzy0iklSNEjR8DdgG3GBm55vZeQSjKZ4Cro52MrPVZjZuZh+Ltrn7/QTDLb9oZu82s9cQDLdcC/xNFa8hlvkdLXztovWsXdLFB757P3uGqhPTDGdytDQZrc3l/5OIujyGqzRXw8+37OHzP32C3z51BR970zpSKWUWRETKoSGCBncfBs4BngC+DXwH2Aqc4+5Debsa0MTh1/Uu4OvAFcCPgJXA6939vgo3fVbmd7Tw5be9mMGxcf6/639dleccGRuvSD0DHCyErEamYSQzzl/8+4OsXdzFJ3/nBQoYRETKqFEKIXH3HcAF0+yzDQ4vcHf3NPDn4a0hnLC8hw+99ng+c9Pj3L15D2cet6SizzecyVVkYic4OGFUNYKGq+98kmf6R7nufS+rWBAkInKkaohMw5HqD85cy9ELOvi7/3qUiYnKTiw0PDZekTkaIK97osJBw+7BUa6+awtvfMHzWL/msOk3RERkjhQ01LH2liY+9NrjefjpAe7c1FfR5xrO5Ois0DfzKBipdE3Dt36+nbHxCS4998SKPo+IyJFKQUOdO//Uo1nW08Y37t5W0ecZGRufHBpZbm3NKZpTVtFMw0hmnGt/sZ3XrVvO2iVxFj4VEZGZUtBQ51qbU1x4xmrufKKPLX1D0x8wS8OZXEWGW0IwF0Wll8e+8YFdHBjJ8u5X9k6/s4iIzIqChgbw1tNX0ZQyrr9vZ8WeY3hsvCITO0W625onJ5CqhOvve5pjl3axfvV0S5GIiMhsKWhoAEt72jjzuCXc8MAuvEIrLY5kxitW0wDQ2dY8OYFUuT21b4RfbtvH/37xMZrtUUSkghQ0NIjzX7SCnfvT3LdjuqU2Zmd4LFfRTENXWzNDY5XpnrjhgacBOP/UFRU5v4iIBBQ0NIjXnbyc1uYUP3ro2bKfOzfhpLO5is5r0NXaxEiFuidufuQ5Tlu1gGMWFi4vIiIi5aSgoUH0tLfwst7F3P747rKfO50NV7isUCEkBOtPVKKmYffAKA/t7Oe1Jy0v+7lFRORQChoayDnPX8bWPcNs3TNc1vNGQyE7KzTkEoKVLisxeuK2x4Ig6jUnLSv7uUVE5FAKGhrIq08MPhhvf6y82YbhyWWxG68Q8tbHdnP0gg5OXN5T9nOLiMihFDQ0kFWLO+ld2sUdT5R3dsgoA1DJmoZKDLnMTTj3PrmXV52wRKMmRESqQEFDg3n5sYv51bZ9jOcmynbOg5mGynVPdLY2MZqdIFfGNTQe2TXA4Og4Z/QuLts5RURkagoaGsxL1y5mOJNj466Bsp0zyjRUcp6G7mjRqjJ2Udz75F4AXqagQUSkKhQ0NJiX9garN/5i696ynXOwKpmGIGgYKeNcDfc8uZfepV0sm9detnOKiMjUFDQ0mGU97fQu6eIXT+4r2zmHRoOgoae9pWznLNQVjswoV11DbsL5n637eOlaZRlERKpFQUMDOn3tIv5n276yTSk9NJYFoKe9kpM7hZmGMnVPbOkbYnBsXGtNiIhUkYKGBnTaqgUMjI6Xbb6GodFxzIJixUqJRmZEWY25uj+cTvvUVQvKcj4REZmegoYG9KKVwQflgzsPlOV8g2PjdLc1V3TYYpTFGCxT98QDTx1gfkcLaxd3leV8IiIyPQUNDej4ZT10tjbx4FP9ZTnf0Og4PRUcOQEHR0+UL9NwgBetXEAqpfkZRESqRUFDA2pKGaccPZ8HnipTpmF0nO4K1jNAXqZhNDvncw2PjfPEc4OculJdEyIi1aSgoUGdtnIBj+waYGx87kMYh8LuiUqKgpJyjJ7YuGuACYdTV86f87lERCQ+BQ0N6pSj55PJTbDpuaE5n2twbJzuCg63BGhrbqK1OVWWmoZHdgXdMievUNAgIlJNChoa1LoV8wB49Jm5zww5NJqteE0DwLz2ZgbLUNOwcdcAi7taWdbTVoZWiYhIXAoaGtSaxV10tDTxSDmChip0T0C4aFWZgoZ1K+ZpkSoRkSpT0NCgmlLG85/XwyNlWINiaHS8ohM7RXraW+ZcCJkZn2DT7sHJTIuIiFSPgoYGtu5583j0mYE5zQyZm3CGM7mKj56A8iyPvWn3INmcq55BRKQGFDQ0sJOeN4+B0XGePpCe9TmiD/FqdE/0lKGm4bFnBgFY97yecjRJRERmoCGCBjNLmdllZrbNzEbN7EEzuyDGcfPM7GNm9nMz22tmB8Kff7sa7a60k54XFUMOzvocUdBQje6J7jIEDVv6hmhOGas1E6SISNU1RNAAfAL4OPBl4LeAe4Hvm9kbpjluFfDHwJ3AhcBbgCeAH5jZ+yvW2io5fnk3AJt3z37YZVSY2N1W2SGXAPPKUNOwpW+I1Ys7aWlqlD9dEZHkqPzXyzkys2XApcCV7v65cPPtZnYccCXw4xKHbwV63X0kb9tPzGwl8GHgK5Voc7XMa29h+by2uQUN4QqX1axpcPdZj3zY0jfMsUu7y9wyERGJoxG+rp0LtALXFmy/FniBma2d6kB3Hy4IGCIbgBXla2LtHL+sh827Z989MTha3ZqGCYeRzOxmsczmJti+d5hjlyloEBGphUYIGk4GxoDNBds3hvfrZnHOVwGPzaVR9eK4Zd1s3j006xEUUdBQrZoGmP1U0k/tGyGbc3qXqJ5BRKQWGiFoWAQc8MM/FfflPR6bmV0CnAH8Xal9zGyDmW3o6+ubUWOr7bhl3QxncjzTPzqr46s7eiKom5htXcOWvmEAZRpERGqk6kGDmb3WzDzG7Y4KPPfZwJeAb7n7d6baz92vcff17r5+6dKl5W5GWR23bG7FkENVzDREU1XPdgTFlr7gGo9doqBBRKQWalEI+XPgpBj7RbUI+4EFZmYF2YYow7CPGMzsJcCNwG3Au2O2te4dHwYNm3YP8aoTZh7gRAtIdbVWp6YBZh80PNk3xJLuNuZ3Vn6kh4iIHK7qQUNYmDiTeoKNQBtwLIfWNUS1DI9MdwIzewHwE+AB4AJ3n9u4vzqyuLuNhZ0ts840DKSz9LQ3k0pVfh2HudY0BCMnVM8gIlIrjVDTcBOQBd5esP1C4GF331rqYDM7Hvgp8CTwJnef/fSJdWouIygOjGSY31Gdb+5zqWlwdzbvHlI9g4hIDdX9PA3uvtvMrgIuM7NB4D6CSZrOAc7L39fMbgVWu/tx4e/LCAKGVuBvgHUF8wPc7+5jlb+Kyjp2WTf/9fAzs5r/oD+drWLQMPvuiX3DGfrTWc3RICJSQ3UfNIQ+CgwBHwSOAh4H3uzuPyzYr4lDr2kdsDr8uXBfgLXAtrK2tAaOX9bNd0ey7B3OsKS7bUbH9qezLKhSjUB3azMpC55zpiZHTqh7QkSkZhoiaHD3HHBFeCu139kFv98BVL6zvsbyR1DMJmg4an57JZp1mFTKmN/RwoGR2QQN4cgJZRpERGqmEWoaZBrRGhSbZlEM2Z8er1r3BMCCzlYOzCbTsHuItuYUKxZ0VKBVIiISh4KGBDhqXjvdbc1smWHQ4O70pzPMq2LQEGQaMjM+7sk9w6xd0kVTFUZ5iIhIcQoaEsDMOHZp14yHXaazObI5Z0FHa4VadrgFnS2zrGnQyAkRkVpT0JAQvUu72bpneEbHRB/e1eyemN8x86BhNJvjqX0jqmcQEakxBQ0JsXZJF08fSJOewQqStQgaFsyiEHL73hEmXCMnRERqTUFDQvSGH6gzyTb0j9Qg09DZysBoltxE/FU5NXJCRKQ+KGhIiN5wEaeZBA3RKIZqzdMAQabBfWazQkYFnr3KNIiI1JSChoRYs6QTCBZ1iqsm3RNhgDKTLootfUOsmN9OZxUW1RIRkakpaEiIztZmVsxv58kZZBoGwqChmkMuJ4OGGRRDPrlnWCMnRETqgIKGBOld2j2joKE/ncUMetqq9w1+fji8M+5cDe7Olt1DqmcQEakDChoSZO2SLp7sG8I9XpFhtFhVNZbFjkRdIXGHXT43MMZwJqeREyIidUBBQ4L0Lu1icHScPUPxvsUfGKneCpeRqHsibtCgkRMiIvVDQUOC9C6d2QiKA+ksC6ocNERBStxCyMmgQTUNIiI1p6AhQXqXBCn8uCMo9g6NsXiGq2LOVUtTiu625vhBw+4hutuaWdZT3XaKiMjhFDQkyIoFHbQ2p2IXQ+4ZGmNJd/XWnYjMZNGqLX3DHLu0CzMtVCUiUmsKGhKkKWWsWdzJk33TBw3uzt6hTNUzDQBLulvZMxwvaHiyb2iy20VERGpLQUPC9C7p5sk903dP9KezjE84S2oSNLTRNzg27X7DY+Ps6h/VyAkRkTqhoCFhepd2sWPvCNncRMn9ohEWteieWNrTxp6h6YOGqKBTIydEROqDgoaEWbuki/EJZ+f+dMn99oYf2ou7apNp2Ds0Nu2iVRo5ISJSXxQ0JEzU/z/dCIrJTENPbTINEw77pymG3LJ7iJTB6sWdVWqZiIiUoqAhYaL+/+mKIfcO1zbTAExb17Clb5hVizppa26qRrNERGQaChoSZkFnKws7W6YddrlncAwzWNRVm0wDMG1dw5Y+rTkhIlJPFDQkUO/S7um7J4YzLOpspamK605EoqChVKYhN+Fs1eqWIiJ1RUFDAvUu6YqVaajFcEs4OGKjVKZh14E0Y+MTk7NciohI7SloSKC1S7voGxxjcHTqqZr3DmdYXIPhlgDdbc20t6RKZho2h5kSTewkIlI/FDQkUO+SaATF1NmGYArp2mQazIwl3W0lV+Pc9NwgACcsV9AgIlIvFDQk0PHhB+2m3cXrGtydZ/tHWT6vdotALe0pPSvk488OsaynjQWdtcmGiIjI4RQ0JNDqRZ20Nqd4/NmBoo/3DY0xNj7BykW1m/9gaXcbuwdHp3x80+5BTljeU8UWiYjIdBoiaDCzlJldZmbbzGzUzB40swtmcZ5eMxsxMzez4yrR1nrQ3JTi+GXdPP5c8UxDNFvkMQs7qtmsQxy9sIOd+9O4Hz4r5MSEs+m5IQUNIiJ1piGCBuATwMeBLwO/BdwLfN/M3jDD8/wD0F/eptWnE5f3TJlpeGrfCADHLKxdpmH1ok5GMrmidQ0796dJZ3OqZxARqTN1HzSY2TLgUuBKd/+cu9/u7u8FbgeunMF53gacBny6Mi2tLyce1cNzA2McKDJVcz1kGlaFU0PvCAOYfE9ERZBHKdMgIlJP6j5oAM4FWoFrC7ZfC7zAzNZOdwIzWwhcRRB8HCh7C+vQieEH7uPPDh722M79aRZ3tdLZ2lztZk1aFdZTPFUkaHg8DBqO18ROIiJ1pRGChpOBMWBzwfaN4f26GOf4DPCYu3+7nA2rZ1HQ8FjRoGGkplkGONg1sn3v4UHDr3f2s3pxJz3tLdVuloiIlFC7r5rxLQIO+OEVc/vyHp+Smb0SeCdB10QsZnYJcAnAqlWr4re0jhw1r50l3a08uPPwxMrO/WnWPW9eDVp1UHtLE0fNay/aPfHQzgP8xpqSL6uIiNRA1TMNZvbacPTCdLc7yvBcrcDVwBfc/ZG4x7n7Ne6+3t3XL126dK7NqAkz49SVC3jgqUODhokJ5+n9aY5ZVNtMAwR1DTv2HToBVd/gGLv6R3nRMfNr1CoREZlKLTINPwdOirFf9BV0P7DAzKwg2xB9Fd3H1D4ELAS+ZGYLwm3RkIEeM+tx98Pz9wlx2qqF3PLobvpHsszvDFL9Tx9Ik8lNsLKGIyciqxZ18t+b+g7Z9lCYGXnhMQuKHSIiIjVU9aDB3UeAx2ZwyEagDTiWQ+saolqGUhmEdcBRwNNFHrsPeBA4dQZtaSinrgw+eB/YeYCzTggyJg9OfijX/pv82iVdXPernQyMZpkX1i88uLOflMEpR9e2+0RERA7XCIWQNwFZ4O0F2y8EHnb3rSWOvRJ4dcEtGnJ5IfDu8ja1vrzwmPmYwQM7DnZRPLDjAG3NKZ5/VO0/lKOg5r7t+ye3/Wr7Pk5Y3lPTkR0iIlJc3QcN7r6bYLjkZWb252Z2tpl9FTgHuCx/XzO71cw25x37mLvfkX/jYJbjF+6+oUqXURM97S08/6h5/GzzwS6AB546wClHz6e1ufYv/WmrFtCcMn65NehhGhzN8sut+yazIiIiUl9q/8kRz0eBK4APAj8BzgTe7O4/LNivicYYEVI1rz/5KDZs389zA6NkcxP8+un+yW/4tdbZ2swpR8/nf7YFQcNdT+whm3Neu255jVsmIiLFNETQ4O45d7/C3Ve7e5u7v9Ddryuy39nuvmaac33D3c3dC+d9SKQ3vvAo3OGmh5/l4af7GRuf4LRV9RE0AJy+dhEPPtXPaDbHLY8+x8LOFl68amGtmyUiIkU0RNAgs3fcsh6OX9bNd3+5g8/d/Dg9bc28/NgltW7WpJf1LiaTm+BTP36U/3r4GV5z0nKaUlbrZomISBEKGo4Al557Ipt3D3H35r38xetOYFFXa62bNOmsE5Zy7snL+dY925nX3sJfnXtirZskIiJTUP//EeDck4/i6+96CXc90ceFZ6yudXMOkUoZn3/zqSz58aO85SUrWTavvdZNEhGRKdjhszNLvvXr1/uGDYkeZCEiIjLJzH7l7uuLPabuCREREYlFQYOIiIjEoqBBREREYlHQICIiIrEoaBAREZFYFDSIiIhILAoaREREJBYFDSIiIhKLggYRERGJRUGDiIiIxKKgQURERGJR0CAiIiKxKGgQERGRWLTK5TTMrA/YXsZTLgH2lPF8taRrqT9JuQ5IzrUk5TpA11KPKnEdq919abEHFDRUmZltmGrJ0Uaja6k/SbkOSM61JOU6QNdSj6p9HeqeEBERkVgUNIiIiEgsChqq75paN6CMdC31JynXAcm5lqRcB+ha6lFVr0M1DSIiIhKLMg0iIiISi4KGKjCzlWZ2nZn1m9mAmV1vZqtq3a5SzOxsM/MitwMF+y00s38ysz1mNmxmt5jZC2rY7mPM7O/N7B4zGwnbvKbIfu1m9lkze8bM0uH+ryqyX8rMLjOzbWY2amYPmtkFdXYtxV4nN7NT6+FazOx3zew/zGx7+G/9uJn9nZn1FOwX628p7mtXi+swszUlXo8F9XAd4XOfa2a3mdmzZjZmZjvN7N/NbF3BfrHeu2r5PhDnWuK+n9X6Woq05aawnVfMpo0V+Rtzd90qeAM6gU3Aw8BvA+cDvwa2AF21bl+Jdp8NOPAB4Iy82/q8fQz4GbATeCvweuBOgjHDx9Sw3c8BPwZ+El7DmiL7fQc4ALwHeA1wPZAGTi3Y75PAGHAp8GrgamACeEMdXYsDXy94nc4AOuvhWoB7gX8H3g6cBXwo/Le/F0jN9G8p7mtXo+tYE74enyryejTVw3WEz/1W4LPA74bX8g5gIzBAMEYfYr531fp9IOa1nM0072f1cC1FruuZsN1XzKaNlfgbq9o/wJF6Az4I5IDj8ratBcaBP691+0q0O/pP9toS+5wf7vPqvG3zgX3Al2rU7lTez++myAct8KJw+7vytjUDjwM35m1bRvAh+7cFx98KPFQP1xI+dsibyhTnqtm1AEuLbHtn2O5zZvK3FPe1q+F1rAl/f/c056rZdZRo04lhm/4i/D3We1edvg8UXsu072f1dC3AQuBZgqCgMGio6f8VdU9U3nnAve6+Odrg7luBuwle/EZ2HrDL3W+PNrh7P/Cf1Oja3H0ixm7nAVnge3nHjQP/BpxrZm3h5nOBVuDaguOvBV5gZmvn3uKpxbyWuGp2Le7eV2Tz/4T3R4f3cf+W4r52ZRfzOuKq2XWUsDe8Hw/v47531d37AIdfS1z1ci2fBh529+8Weaym/1cUNFTeyQTpvUIbgXVFtteb75hZzsz2mtm/FvRnlrq2VWbWXZ0mztjJwFZ3HynYvpHgg/W4vP3GgM1F9oP6ev3+KOzPHQn7d19Z8Hi9XctZ4f2j4X3cv6W4r121FF5H5O/MbDysBbixSH9zXVyHmTWZWauZHU/QXfUsEH1QxX3vqov3gWmuJVLq/Qzq4FrM7BUEGaz3T7FLTf+vKGiovEXA/iLb9xGkoOpVP/B5grT4OcAngNcC95jZsnCfUtcG9Xt907V7Ud79AQ/zeiX2q7VrgT8meH0uARYDt5nZ2Xn71M21mNnRwOXALe6+Ie/54/wtxX3tKm6K6xgj+MB6L0HdyKXAC4Cfm9lJeYfXy3X8gqDNTwAvJOhm2Z3XhjjvXfXyPlDqWuK8n0GNr8XMWgn+fj7n7o9PsVtN/680z+YgST53vx+4P2/TnWZ2F/BL4E+Bv65Jw+Qw7v6OvF//28xuIPgmcgXwitq0qrjwW9ANBGnjd9W4ObM21XW4+zPA+/J2/W8zu4ng291HgQur2c4Y3gHMA3oJApyfmtkr3H1bTVs1O1NeSwO9n/0V0EFQtFyXlGmovP0Uj06nigLrlrvfRxDFvyTcVOraosfr0XTt3pe33wIzs2n2qyvuPgj8iIOvE9TBtZhZB0G/ay9wrrvvLGhfnL+luK9dxUxzHYdx96cIqt0LX4+aXgeAuz/q7r8I+85fA3QDHwkfjvveVRfvA9NcS7H9C9/PoIbXEnaVfBT4P0CbmS3IG6Yb/d40gzZW5G9MQUPlbSToWyq0Dnikym0plyjFXeradrj7UPWaNCMbgbVm1lmwfR2Q4WC//0agDTi2yH5Q/69ffldETa/FzFqA64D1BEM8f12wS9y/pbivXUXEuI5SCl+Pml1HMe5+IHzeqK877ntX3b0PFLmWkrvn/VzLa+kF2gm6G/fn3SDInOwn6Oqq6f8VBQ2VdyNwhpn1RhssmKDnzPCxhmFm6wmGMv0y3HQjcLSZnZW3zzzgf1Hf1/afQAvwe9EGM2sG3gLc7O5j4eabCKqP315w/IUElc1bq9DWGQtfgzdx8HWCGl6LmaUIxoufA/y2u99bZLe4f0txX7uyi3kdxY5bRdBNlP961Ow6pmJmy4HnE8zDAPHfu+rufaDItRTbp/D9DGp7LQ8Q1MEU3iAIJF5N8EFf2/8r1Rp3eqTegK7whf41wXCY84AHgSeB7lq3r0S7v0PQJ/6/Cd4k/4Jg8pAdwJJwnxTwc+Ap4PcJhvXdQZD2WlnDtv9uePsqwbeIPwp/Pytvn38jiNzfTZDKvA4YBV5ccK4rw+1/TjDW+6sEEyK9qR6uheAbyNeAt4Xtuyj8W8sAr6yHa8lr+xUcPuHRMTP9W4r72tXoOj4PfAF4M8Gb/PuA7QQT7JxYD9cRPvcPCNLg54ftfC/wWNjOE8J9Yr13zeS1q+G1TPt+Vg/XMsX1Fc7TUNP/K1X/BzgSb8Aq4D8IZigbBP4vRSbpqacbcBnwEEHVcTb8A70GeF7BfouAfwn/YEcIJgt6UY3b7lPc7sjbpwO4imBY1ihB5fXZRc7VRFAktZ2gMvsh4Hfr5VoIvl3cHb4BZgnGp98InF4v1wJsK3EdH5/p31Lc164W1wH8AcHcDfvD1+NZ4F8pCBhqeR3hc38Y+BXBB+sIwYQ/Vxe+L8V976rl+0CcayHm+1mtr2WK6zskaJhJGyvxN6ZVLkVERCQW1TSIiIhILAoaREREJBYFDSIiIhKLggYRERGJRUGDiIiIxKKgQURERGJR0CAiIiKxKGgQOUKZmce4bTOzNeHPF9e6zREzO9rMhsOpgKvxfGZm95vZX1Xj+UTqlSZ3EjlCmdkZBZt+QDBN8Mfzto0RLE50GrDF3fuq07rSzOxfgGXu/qYqPufvEMzCd6y71+UKpyKVpqBBRAAws23Az9z9wlq3pZRwMaKngN9x9x9V8XmbgJ3AF9z9M9V6XpF6ou4JESmpWPeEmX3DzHaa2Xoz+7mZpc3scTN7Y/j4n4ddGwNmdoOZLS04Z7OZXWZmj5nZmJntMrPPm1l7jCZdTLAOwk8KznmHmf3MzF5vZg+EbbrfzF4aPt+nzOwZM9sXtr+roD2fMLMtZjZqZnvCc70i2sfdc8D3CRb/ETkiNde6ASLSsOYB3wI+B+wCPgr8h5l9BTgBeD+wHPgi8BWClR8j1xIstvVpghX7TgI+AawBLpjmeV8P3OPu40UeOw74LPBJYAj4DMECXjcSvN9dHD7XZ4HdQFSj8GHgz8JreCC8tvUECwPluwv4gJn1uvuT07RTJHEUNIjIbPUA73P3uwDMbBdBTcSbgHXhN3PM7BSCD9omd8+Z2SuBtwAXufu3wnPdYmb7gGvN7FR3f6DYE5qZAS8lWH66mMXAy6MPdDNLATcAa939teE+PzGzVwG/x8Gg4WXAze7+/+ed6z+LnP/+8P4MgiWiRY4o6p4QkdkajgKG0GPh/S1RwJC3vRl4Xvj764EMcF3YLdBsZs3AzeHjryrxnAsIlvudqiDziYIMQNSmnxTs9xhwTBiEQLCc9RvM7JNm9goza53i/NHzrijRRpHEUtAgIrN1IP8Xd8+EP+4v2C/aHtUrLANagWEgm3fbHT6+uMRzRucYm+LxqZ672PZmoCn8/VPA3wDnAf8N7DWzr5vZkoLj0uF9R4k2iiSWuidEpNr2AqPAK6d4fNc0xwIsLGeD3D1LUF/xaTM7iqCL5Sqgk6ArJRLVOOwp5/OLNAoFDSJSbTcRFB7Od/dbZ3Kgu2fMbCvQW5GWBc/xLPBPZvYG4JSCh9eG949X6vlF6pmCBhGpKne/w8y+S1DTcBXwS2CCYOTEG4APu/sTJU5xF3B6OdtkZjcQFHHesz8b8wAAANtJREFUR9CVcRpB7cXVBbu+lKAr5d5yPr9Io1DQICK1cCHwAeAPCIY5jgHbCAoWn5vm2O8B7zSzNe6+rUztuYtgNMX7CbokdhAM1/xkwX5vAm5095EyPa9IQ9GMkCLSUMJhlJuAr7v7FVV83hUEM1G+bqbdKiJJoaBBRBqOmb2doFBxbbW+9ZvZF4AXufs51Xg+kXqk7gkRaUT/ChxNUAfxSKWfLJzP4Vngmko/l0g9U6ZBREREYtHkTiIiIhKLggYRERGJRUGDiIiIxKKgQURERGJR0CAiIiKx/D/ndax0DnIHXAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "src.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our PMLs, we will need some damping parameter. In this case, we will use a quadratic taper over the absorbing regions on the left and right sides of the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Damping parameterisation\n",
    "d_l = (1-0.1*x)**2  # Left side\n",
    "d_r = (1-0.1*(grid.shape[0]-1-x))**2  # Right side\n",
    "d_b = (1-0.1*(grid.shape[1]-1-y))**2  # Base edge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now for our main domain equations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from devito import Eq, grad, div\n",
    "\n",
    "eq_v = Eq(v.forward, v - dt*grad(p)/density,\n",
    "          subdomain=grid.subdomains['main'])\n",
    "\n",
    "eq_p = Eq(p.forward, p - dt*velocity**2*density*div(v.forward),\n",
    "           subdomain=grid.subdomains['main'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also need to set up `p_i` to calculate the integral of `p` over time for out PMLs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "eq_p_i = Eq(p_i.forward, p_i + dt*(p.forward+p)/2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And add the equations for our damped region:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Left side\n",
    "eq_v_damp_left = Eq(v.forward,\n",
    "                    (1-d_l)*v - dt*grad(p)/density,\n",
    "                    subdomain=grid.subdomains['left'])\n",
    "\n",
    "eq_p_damp_left = Eq(p.forward,\n",
    "                    (1-gamma*velocity**2*dt-d_l*dt)*p\n",
    "                    - d_l*gamma*velocity**2*p_i\n",
    "                    - dt*velocity**2*density*div(v.forward),\n",
    "                    subdomain=grid.subdomains['left'])\n",
    "\n",
    "# Right side\n",
    "eq_v_damp_right = Eq(v.forward,\n",
    "                     (1-d_r)*v - dt*grad(p)/density,\n",
    "                     subdomain=grid.subdomains['right'])\n",
    "\n",
    "eq_p_damp_right = Eq(p.forward,\n",
    "                     (1-gamma*velocity**2*dt-d_r*dt)*p\n",
    "                     - d_r*gamma*velocity**2*p_i\n",
    "                     - dt*velocity**2*density*div(v.forward),\n",
    "                     subdomain=grid.subdomains['right'])\n",
    "\n",
    "# Base edge\n",
    "eq_v_damp_base = Eq(v.forward,\n",
    "                    (1-d_b)*v - dt*grad(p)/density,\n",
    "                    subdomain=grid.subdomains['base'])\n",
    "\n",
    "eq_p_damp_base = Eq(p.forward,\n",
    "                    (1-gamma*velocity**2*dt-d_b*dt)*p\n",
    "                    - d_b*gamma*velocity**2*p_i\n",
    "                    - dt*velocity**2*density*div(v.forward),\n",
    "                    subdomain=grid.subdomains['base'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add our free surface boundary conditions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def freesurface_top(p_func, v_func):\n",
    "    time = p_func.grid.stepping_dim\n",
    "    pos = int(max(p_func.space_order, v_func.space_order)/2)\n",
    "    \n",
    "    bc_p = [Eq(p[time+1, x, pos], 0.)]\n",
    "    bc_v = [Eq(v[1][time+1, x, i], -v[1][time+1, x, 2*pos-i]) for i in range(pos)]\n",
    "    \n",
    "    return bc_p + bc_v\n",
    "\n",
    "bc = freesurface_top(p, v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And our source terms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "src_term = src.inject(field=p.forward, expr=src)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct our operator and run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "tags": [
     "nbval-ignore-output"
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Kernel` run in 0.04 s\n"
     ]
    }
   ],
   "source": [
    "from devito import Operator\n",
    "\n",
    "op = Operator([eq_v, eq_v_damp_left, eq_v_damp_base, eq_v_damp_right,\n",
    "               eq_p, eq_p_damp_left, eq_p_damp_base, eq_p_damp_right,\n",
    "               eq_p_i]\n",
    "               + src_term + bc)\n",
    "\n",
    "op(time=time_range.num-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is important to remember that the ordering of equations when an `Operator` is created dictates the order of loops within the generated c code. As such, the `v` equations all need to be placed before the `p` ones otherwise the operator will attempt to use the updated `v` fields before they have been updated.\n",
    "\n",
    "Now let's plot the wavefield."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "scrolled": true,
    "tags": [
     "nbval-skip"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAGDCAYAAAAyM4nNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8GearUAAAgAElEQVR4nOydfZgcVZX/PycMw5CEmARCACEkGAWJyIsDi4ousr6gIoggoKAoEsQVFBVWWFAR8AcuiqgsQvB1wVUWFAFFQCQgIhGjgBIEjSTyYtDEhJcQ8kbO74/q7jp1p29NTXXPTPXM+TzPPHOr6t6q29VVfaq+595zRFVxHMdxHKfzGTPcHXAcx3Ecpz24UXccx3GcEYIbdcdxHMcZIbhRdxzHcZwRght1x3EcxxkhuFF3HMdxnBHCqDLqIrKdiFwtIk+JyNMi8kMRmTbc/XIcx3GGFhHZVkS+KiJ3icgqEVERmV6w7RgROU1EFovIahG5T0QOidSdLSIPisgaEXlIRI5v5+cIGTVGXUTGArcCOwFHA+8BXgzMFZFxw9k3x3EcZ8iZCRwGrADuGGDbs4EzgYuANwPzgKtE5C22kojMBi4FfgDsD1wFXCwiH2qp5znIaAk+IyIfBS4AdlTVhbV1M4A/A/+hqhcMZ/8cx3GcoUNExqjqhlr5WOAyYIaqLu6n3ZbAo8B5qvoZs/7nwBRVfXltuQv4G/BTVT3a1PsmcCCwtaqua++nGkVv6iQncV7doAOo6iLgTuCgYeuV4ziOM+TUDXoJ3gR0A1cE668Adqm9LAK8EpjSpN7lwObAPiWPn8toMuqzgPubrF8A7DzEfXEcx3E6k1nAGmBhsH5B7f/Oph70tTthvbYymoz6ZBLfSchyYNIQ98VxHMfpTCYDT2pf3/Vys93+D+1OWK+tdA3GTkcKInIccBzAuHHjXrHTTjsNc48cx3Gqw+LFi1m2bJm0e78zRXRVi/tYkrwRrzar5qjqnBZ3W3lGk1FfQfM38tgbPLULYA5Ab2+vzr/77sHrneM4TofRu9deg7Lf54BWh4d/Glaram87+hOwApgoIhK8rdffvJebepDYnSU59drKaJLfF5D6OCw7Aw8McV8cx3GczmQBsAnwomB93Uf+gKkHfe1OWK+tjCajfh2wt4jsUF9RCzTw6to2x3EcpyKMafFvELkRWAccGaw/Cri/NqsK4C5gWaTecpKZV21nNMnvlwEnANeKyBmAkgQQeJQkOIDjOI5TAYSheeMUkUNrxVfU/r9ZRJYCS1X19lqd9cB3VPUDAKr6DxG5ADhNRJ4BfgccDuxHMnWaWr11IvIpkmAzjwO31OocA5yoqmsH4zONGqOuqs+KyH7Al0jmCQrwc+AkVV05rJ1zHMdxMgyRjHxVsHxx7f/twL618ka1P8vpwErgo8BWwEPAYar6Y1tJVS8REQU+AZwCPAKcoKoXM0iMGqMOoKqPAE3j8zqO4zijC1Xtd+R+szqq+jxwTu2vv/aXMoRq8Kgy6o7jOE5nMJoGfLUTN+qO4zhOpRgqn/pIxI264ziOUzncqJfDjbrjOI5TKfxNvTx+3hzHcRxnhOBv6iXYYJ6FxlAwe9/69Wl59erstp6eRnHV+u5G+ckn403s7rrMtzhuXFqeMiXbZsyyf6QLK80sPruziROzjbbYolFc/mT6ucP+1Bk/PrtsPhrdXeZcxT5AwAZ/7nQCovecvaZsOcRclE+vTK+vJ57IVlthgkdvvHFa3mqrtLzNRBOhfGGQtGvZsrRsr3FzTzF9eqbJ358Z2yg/9VS6Pna/Qfaes7fv5PFmGvTKYNauvTFt2ZB37xX+3WsBv/PL4UbdcRzHqRxtzxIzSnCj7jiO41QKoW+0F6cYbtRLYKWnUKKKylJWPwulMCMVrjYpdq0caJW8cHdWftt227Q8dUrQF7tDW163Li1vvXW2jdn5ypWpNBhTFkOZ0CqN3YE03yBHKh1jdu5S/OgkV+q1F5wth9eUvY560uvY3gYPPphtYrdFFGq26TL38h/+kN340EPNG+24Y1oO/FXPrp/WKFs1P7z/LfYee9nL0vJkTN/CHVid3t7AOa6wocbv9nL4eXMcx3GcEUJ1Hsscx3EcB5/S1gpu1AeT2Ij3HC1t4stS+d02CUfm5g3ubRDK/HY4/dKladkO7c0ZZhv7OHl9saPxx4+3swYiO4ZKSYBOBYmNcs+7KDNupHS1vRUfeyzbxN5z9joOBqyn2PsLYPHi5v2x1/eLX5xpMm77VH5/5pl43yyx0e8sNMfPG/1eUdyol8N/PR3HcZzK4Ua9HH7eHMdxHGeE4G/qjuM4TqVwn3p53Ki3mUy0uZhvOPSpL1mStjF+wZ122qtRDt11dtn61F7wgkglyPr1Nt00LW+ySfOdQcb/Zz/ORmYSaV4grzVr0rJ1eY7tMTsrNEAgfyqhM7LITGMLr48ifvTw3jPX9TLjn7a3Yl7UNnu7ZHZtFzbbLNvIOrjtgax/O3CWT91tt0Z5663T6JI2ulz40WbOTMsTHnsgXbj33rRs572FVHTaaHV60lm4UXccx3EqheAR5criRt1xHMepHB5Rrhxu1AdAXZoqlczATiEJ5cT77kvLv/51ozh531Sae8v++2ea/O3JNCqW3fXkHpNg4tngOFYOtDK7LQcy3YauVAK0Kr1VGu2MuBAradopOl1dqbjWnTeFrWDiF2cEE94vRST3wI309Or0OrazOa0SHl5e9laYOjWya7tgdfCwb8bFliGcanb//Y3ia/dJpfiZM9P7JVTSu395a7rw4x+nZXuegt8P+1tQRnKvkkzvZPFfScdxHKdS+EC58rhRdxzHcSqHG/VyuFEvQVHpKTMS3sp0NiEzZEepf//7zcuvf32myTavfnW6YKNSWW0uHMlus73YbUa/X7s++GyRAfNWybcKYl7e99hI+K7xqTQKMGb9WpxRTt6UithsCXsd90zIbHpicVq2I8ktoaxtr3F7y9p61j01xmZTCRvFsjOFM1TsiPW77moUt7H1FizItpk7t/n+Dj44Ldt7n+D8VNTD5Ua9HBX6Ch3HcRzH5fdW8PPmOI7jOCMEf1MfRKxKuHJl+vw08WUvz9Qbs+uujfLqyy5rlO81w8W7v/e9TJsdzPLEPfdMNxx6aFoOJPsNM1/SKNuYF1YKD/M8RFT6Psp+nVC+iw1Ufu65eJsygWksRZtUSWp0cgLO5H2h9qI0evmyIAGSVb+tu8gGUZo0KdvGKtZWSbfXfmZ2x3NZyf95I3FP3G2HRnnsskfSSnfckT3ojTem5WuuaRQfNgf6W7ZFZi7ObqY8obe3UV610x6ZNvb8WHdCle4Jf+MsR4W+QsdxHMdx+b0V3Kg7juM4lcMjypXDH4Ycx3EcZ4RQ2Td1ETkUeBfQC2wJPAL8EPh/qvpMrc50YFFkF5NUtTG/Q0R6gLOBo4CJwL3AJ1X1F4P0ETJY39tvfpPdtvvsDzfKPT//eaO83vjUfhfs7/emvJvZ4W42XJadkwOMMRGvVq9O/X0rVqR17LS1EOtviwXvCt2fXRH3uJ3eFkakKxJtLozqVybCVcxVWyW/YqdS6twW9aNbjE99+ZPpNfBE4FO3M73WrUvLNjJiONPULk+e2DyKpB0r889/xo9p/dYv2cLcMIsXZxuZsTK3mNX2/g8nfFpv+YTXvKZR3nDSxxvl24yrHmD69OZ9qxJDESZWRLYDvgS8gUQcuAU4SVUf6afdmcBnIpvXqGqPqbsY2L5JvYNV9Uclup1LlX++TiYx5P8JPAbsDpwJvE5EXqWq9i47F7guaP9MsPwN4K3AKcDDwIeBm0Tklap6L47jOE4lGAqfuoiMBW4F1gBHAwqcA8wVkZer6rM5zb8OBI9KjKutC20RwE0k9svyUIlu90uVjfrbVNW8dnK7iCwHvgPsS/Jl1HlYVefFdiQiuwLvBo5R1W/V1t0OLADOAg5sc98dx3GcFhgC3/BsYAdgR1VdCCAivwf+DHwQuCDWUFUfI3nZbCAi7yGxqd9p0mRZno1qJ5U16oFBr1PXmV84wN0dCKwDrjT7Xy8i3wdOFZFNVHVNtHUORWVgqyb+5S/ZbT/5SVo+9sIfNsr7LHtto7w6mPryoClbKb7byHk7X3RR9kBG+3zJ29/eKP9t/ORGOcwvYbHT0GIyfSiv2mXrgshElMtpQ8+YaL3BosQsugwjQb5v9RwUpdXkSKtIExvlBWqzn8e6e6yHKvBWZZfz5n1GVlv5fofp5nP+6La0fOWVWKxkaO9xk6aJl5DlLeZAa29JvYlnnJrWCQLKZeT32PWa990MdkKXIRr9fiAwr27QAVR1kYjcCRxEjlGPcDTwd5K38mGj0wbK/Wvt/x+D9eeKyHoReUpErhORXYLts4BFqroqWL8A6AaC9EqO4zjOCGcWcH+T9QuAnQeyo5pv/nXAd1W12SPx20RklYisEZF5IvL2JnXaQse8U4jIC0mk8ltUdX5t9RrgUuBmYCmwE4kP/lcispeq1o3/ZGAFfVlutjc75nHAcQDTpk1rx8dwHMdxCtCGN84tRGS+WZ6jqnPMcp5dmNRkfR5HkXS5mfR+PYnKvAiYCpwAXCMi71HVKwZ4nH7pCKMuIuOBa0lSD7y/vl5VlwDHm6p3iMiNJE9ap5Oc6NLULoA5AL29vVqXo/KkJytZ2VHcsQhskE2BfP75afnMM1Mp7ZPv+2amzd4f+ECjbEfG2mhTXX/4Q6bNS046KV0wiSO2OeCAdP1uNiYVrJ24ZaMciw6XR6xemPilCFbeLCNx543MbydDJV2POAoOmd8wPp25scx4NfPk99iId3tfhvdoVH42F+9kc01OHh/4rh40Avq3zVh2k6jpkT//OdPEdnsnU7avjducckqmzXX7/FejfILRHMeNS8snnpjtmnUt2NNb1B1Sym0yQNpg1Jepam//1drCe4F7VPX34QZVzZx9EbkGmEcywLvtRr3y8ruIbErypLMD8KbaAIUoqvoo8EvAxE5lBc2fvOpv6MubbHMcx3GGgbpPvZW/AuTZhWZv8M37KrIXyTNYs7f0Pqjq88BVwLYisnXR4xSl0kZdRDYGriaZq/4WVf1DP00sasoLgBm1KQyWnUmmfS7EcRzHqQzS4l8BFpD41UN2Bh4YQFePJhmI/b8DaFNH+68yMCorv4vIGOC7wH7AAUWnA4jINGAfwE7qvx74LPBOak9TItIFHA7cPNCR72UkqvHj0+enF70oW89KYatXp8FnTj01FSVO5VDbhB/84JhG+R0LU/lt1Sc/2Sg/SJa7TYKYbb761UZ526uuSiu9852ZNt0mKUy3lebNiNtV69Oc0nmyeiy2SN6I+40iESgiA5BzjzmQbXVGwkj2jifQxW18JSu523J4TW2ySVqOJSbq4yKzO4kdyMrnv/51tv3116dNzKyUf5gq4eW1r1348pcbxRtmfqRRnh4MsVp3vp0Wnb547rlnGogmTPVuA84MhZReUa4DviAiO6jqw9AIaPZq4NScdg1EpBs4AvhpZMZWszZ12/OIqj7RX/2BUuWfrP8mMcKfA54Vkb3NtsdU9TER+SKJ2nAXyUC5HYHTSBIXfa5eWVXvEZErgQtrb/+LgA8BM4Ajh+LDOI7jOMUQhiSi3GUkg9auFZEzSN6azwYeJRmAnfRFZHvgL8BZqnpWsI8DSOT6ptK7iLyLZHrcDbX9TiUJfLYHScTUtlNlo/7m2v/Ta3+Wz5JE51lAYpzfB4wH/kkSlOazqhpG63k/iaE/hyRM7H3A/qoaRmB1HMdxhpnB9g2r6rMish9JmNjLSZ4lfk4SJtZqPfVnjGZdOppkTNaPm2yD5AVyS+B8EuP/LDCfxPYMynz2yhp1VZ1eoM43gW/2V69W9zng47U/x3Ecp8IMxYCvWoz3Q/qps5iIm15VD+qn7TwSF/KQUVmjPtLoNmkYdtmlO7PNBHdj7lwbB8dGm7o60+aQQ9LB/Rtv/OlGefHj/9Eo73HsWzJtHv7pTxtlOzLwEZP9YprxtQNs+yMzNMFOfdt330ZxrPG1j52ZjePz96XprWl92DY63fPPUwhbL296WiyJTB4x33mZ6WkjwQ/fzml5hc+HrRhJ1ALw97+nZZu4xU5jC/sfixxno751L7MTQslOSbvfxCiZZ4b3mHvqsWAenZ2mY73W0015m/e/H8vvT0rfUXbd1Q5Q+ZQph6kqbK6QTzRKRxyRru0NJnZ1r7dxuCI3j9OR+DfoOI7jVIohChM7InGj7jiO41QON+rlcKM+VBg9cOzKrEx37LFp1LbFi1Mp7Utferep9cVgh2kO9XXr3tYov/CFsxvl44+/IdPiazfd3CjvYKT0+SbcVhgIefGjjzbK07/2tUZ526uNOyAiywNM3TudtDBpZpqKwqaRXhGEebDT4mwkMCuphtOPrGoYmwYXYhN72H3HFMiiyuRoiSjXjs9poy6uJ3VL2VskzI1ul+3sMnvdvOAF2TY2oYmdmTlm3q/Shblzs41uv71R3PCznzXKD5sqVrAPE0tsacp72DB2N6YZO9990asybb63qw0wlk32EidNMjl79oxG+dhj0xpjnwxcC3ZeX5jJpgL4m3p53Kg7juM4lcONejn8vDmO4zjOCEFU2x6lbkTS29ur8+++u/wOrDYY6olWQ9wpTePw+a+liStOPXVRsMOvmHIsyu2OwXI6Mv7aa1Nx8MA/NOL0sPKMMzItbHYC6zSw8XZtquY+OWwPPjgt779/WjYy/drp2QzRj5lhw7Foc6EUHoswV2Ywb1Ep38r3RfZVdYrI6dYdUoZNNy12fDuQ3F4PEJfcrYocRlCbtt6I5nZGh0mu8vRvfoPFyuz2jl1ryja94x7ZQzL2859vlH84M733Djlkian1ObL8leZYMf8jmS2nnLJro/xf55lx9vNNgrLQX2X9EXnZpvqhd6+9mD9/fsGorMWZJaJFnQ8xdoHfDmFCl8rQQT85juM4zmhgiCLKjUjcqDuO4ziVw33D5XD5vSAty+9WWwzldxvkwg4LN4FcHtgqG5ToeJNF/o47bjZbvmfKRqfsw781Sj09aZ71O+/M1trj6v9slJ8+99xG2cry9ijZsDqwjSlbaX68GRWfkegB9tknLRsd9WlSd0Qoy8eSxRRN6FJmxHuro+RjFB29n0fRgD4xYuetzIh36xoJP5t1YVgp3crvJg9RH7bbLi1vY8eiX3NNtuJ30tDc/zAyu016FGbWsB/VxKjJyOwTTzutUf71Qf8v0/41aT4V1q2zs1duI84EU05nv+y++1sb5QsvzLZ47XgT6doGydl887T80pdmG9moO0WzIzVhsOT3XUT0hy3u4yWjVH73hyHHcRzHGSG4/O44juNUDn/jLIcbdcdxHKdyuFEvh/vUC9KyT90SJH7I+NRtsgjrX99ii2wbE8XtusdSL5+NIrV0qY1OBXCtKdtkEdbR+VYskyalEepM7gr+5aY0rfDqz3ymUQ7z2Fo/pb1J7afZIWizzZQp6cLrX5+WrR8+nLM0fXpaNv7CVWbyXXjarQ/XlsvQ6tS5oaKMT7xMUhzrpt1kk7Q8aVK2TXeXmYJlv4TY/DaAhWYK5223peWrrmoU/2bvKbK+cxtbzSZasX5zgL1M2frOf3do6ju3l+eKFZeSJZaN0/qws0m+Jk06qlH++tfT9e/Y6YF0wUSkA7Jz/Kyv3IbOM1Nlgez8vxYuxMHyqb9cRK9vcR/T3afuOI7jOE4n4/K74ziOUzn8jbMcLr8XpK3ye6j12ilu95pcyTYi1MIgapyVzKzM9r73NYoXXGGjUMEnPmGP+01TtlPiwrlQVpp/oyl/sFH6wQ/SOu9Y/3/Z5iec0Cg+sHRpo/yIqRIq3za+lZXp7fS4LW2CDMieA1u2Mn0oQUYk+6dXpxPzbN73NWuyzctMo+sUYvnpbdnK6pCV1rtXLk8XbEg461ICMAmDePzx5m0CKX2DmZJm92avKXP0PtjryErsYy66KFPvhhkfbpRtziLVi00t45PKxd47xzRKX/ziuEytjx/1j3Th299Oy/eHqZYM9jq2idPttb9V4FxoYRqbZTDl9xv6r5bLdqNUfvc3dcdxHKdSeES58rhRdxzHcSqHy+/lcPm9IG2V30NseDQrO1rJzcrykJUx7ehgO0o+yG2+4X2p7Hfqqen688//p6l1edC520w5FqbMjErnqMyWj30szQ9/wRlGFD355EZx9be+lWnzJ1M2YmQmkUZ4w1vJ3ibZsA6ILULJ0cr0u++elk0kv5hEn+zQnGs7mtgkyNjQY1PfZKX52GDvkCJJVPKSy8RGpYeDnsesNhnB7TVpry870jrMtPLnP6flBQvS8j33NIobAindyuQ2SKDNTR66Z9bSHHsNTAu2TXjzm9OFL3yhUfz89Ts3yqeeau8DyLqofhE5qn2f3DfYlkaEO+WU9Eo855y0RvcV9hhkR/Pb826TrtjrE+Kj3O31GiZtadPUi8GS33cV0Zv7r5bLVi6/O47jOM7wI/ibelncqDuO4ziVw416OVx+L8igyu9We7Wypx0VvyjIp/7QQ2nZjoy3bUJN18pxJmrG0/sf1iifeWa2yZe+tMIsXW3Kt5ny08TZ3pRTCXTKlDTIzXnnZVscs6/JZG03XnZZo2glesgGFLGhSuwZCJPNxCT7iZH1E0KNOybN23IYNMjK9LZsdfHQTVBEKg2/61hkHSvphkFdrLRuryM7Qt24fcIR5ssiZXuUUEq3P9w2lYl1m9jR6gCT7Xm3yYBM5KXrFqayOmSv63vusVmLfmLKfyCO7d1rTTm9dz72sWxknTPOSMuTbzGzQqzEHrow7Hdtc55byf1FL8q2mTEjLcck90GKdDRY8vvuInpri/uY7PK74ziO41QDf1Mvh583x3EcxxkhuPxekEGV3y0xKT6USq08aiU8OwL5r3+Nt7FYmS8YMb/27am8aAYNZ+TMdeuuMy3mBjsPguY0sGJ4EMfdyJubbZbmfbcj9s3g+WRvl3wlXTBBs5f/IZVUA6EzIxHb0dax0dWhrNUdKfdEynltbDk8TpEn7w3BshXj1xYoQ1Yaj5Xz2ts+xGT1ML76S0x5rJXSjz++UXxgWxu4BS65JC0bjwyrV1tZ3eRQAMC4qzJzKuzo9WBUOa9ulDbe+MBG2crq9joce2OQATxPZq8Tume2N+6qF784LVuXQ9jGLudNbxgEBkt+30NEb29xHxNcfnccx3GcauAycjncqDuO4ziVw416OSp73kRkXxHRJn9PBvUmicjXRWSZiDwrIreIyC5N9tcjIueLyBIReU5E7hKR14b1HMdxnNGBiGwnIleLyFMi8rSI/FBEwthFsbbN7JOKyG5BvTEicpqILBaR1SJyn4gcMjifqMI+dRHZl8RB+xHgN2bTelWdX6sjwB3AdOAUYAVwGjAL2E1VHzP7+y5JsvBTgIeBD5PMsXqlqgbh2voyZD51S16WkJi/3frNw+QZdjkWkS48jvW322QRhx7aKF7509Rr+uUvZ5vfdZf1Zd5myr815SUUw2Ze3zPYtm+j9LrXpX3+YJp3hsPfZj3nwNVmit6PTe7rO+5oFJ8259P64CE7VSsWDS0MFGeXN0TK7WZMpBzKdHbZjgWwU/+sf9xO94MgYt9rXpOWDzI5w2fPzrT5yiXpaAI7ZuPRR61P3PrKAUxucWwUOBvxMJsoJTu10v7mvqpRes1rbB34cJrPhcPfbKZt/uhHaXme6WdszApk/d72nrK+8nDZTk+LRCwE4hl3hoDB8qm/QkR/1eI+evrxqYvIWOA+YA1wBqDAOcBY4OWq+mze/kVEgW8Dlwabfq+qq0y9zwEnA6eT/PAdAcwGDlBtOW9NHzpBfv+jqoajXuocSDKaZT9VnQsgIncBi4D/IHkgQER2JYnZeIyqfqu27nZgAXBWbT+O4zhORRgCGXk2yZvCjqq6EEBEfg/8mSQN5QUF9vF4jn1CRLYkMejnqWr9sXWuiMwEzgPabtQrK78X5EDgb3WDDqCqTwHXAwcF9dYBV5p664HvA28SkSCJpOM4jjNc1MPEtvJXgAOBeXWDDqCqi0ikoYOirQbGm0gmt1wRrL8C2EVEZvRt0hqd8Kb+XRHZgkTtvAk4VVXrqZNnAc0SDS8A3isi41V1Za3eIiuJmHrdJHNZFlA18mS1WJYOK82FU19isp+V4sOpN1ZStLLjL3/ZKB5ukkgcflEaqQ7gkS32aJS//e0dG2Uz64xHH70ve8xM8gzrGXk4UgbzvMbcueNMOZ0WdESfqXP/0ihNmvTeRtkE28uUg9l+7DHRTI0yOb4zyUwWBlP67Lm259a6QJ55JttmdRiHrQlhFDqbb95Kt1bSDaVfG7Vs1qy0bJLdPLI+je92c3oJANkZXLfckpYXnWCiIZ7wvWyjzNSzyLSvXKxkbhMLvSZTa88903pHmZxDhx+elqf+NXCvmWuck00iGht5z96HRe83uz5MEhST2Yd4qloVGII3zlnAtU3WLwDeWXAfHxKRU0h8P/OAz6jqHWb7LBJ5P5zbW/+R2JlEWW4bVX5Tfwr4InAssB9wNvB64K6apAGJW29Fk7b1CJaTCtYL3YMAiMhxIjJfROYvXbp04J/AcRzHGS62qP9+1/6OC7bn2YVJTdaHXAH8O4ldOg7YHLi1Nh7MHuNJ7Tt4Ldf2tEJlH/lU9R7gHrPqdhH5BXA3ia/8jKYN29uHOcAcSAbKDfbxHMdxnIRkHHQLqC4bzOAzqvoes3iHiFxLohyfA+wzWMftj8oa9Wao6u9E5E+kQ59X0PyJarLZXv+/fU69MDdF9bESXEymC0fJxmRYKwcuCUaiW7n4739Py1aCtCOAg7zv08xxPv2yVP7+9G37NsoPrN410+b730+XrzCeqEWL7Eh6OyECssk4bCS9eyNlsG6uFSvSyGJXXTXBlG1qEXOe+izb8rtNeetMi623To9jT/vmRu2eFFzRMeU1lpsdYIV5//inGSD+mHFWLfmZHS0O2VkItvzryHpzPSRHMmU7VyB3EHEEK2WHMx3SKIO77566dEw+F444Itti8kIjrdvr9XORZEiQPcH2S7BfnC2/8IXZ9tttl5atNJ83kn0UyuxNEWn9869b11+NPPvR7A0+F1V9RrzeciYAACAASURBVER+AnwgOMZEEZHgbX3QbE+V5fc86idnAYnPImRn4JGaP71eb0ZtCkNYby3xWKaO4zjOcNDV1dpf/+TZjwearC+KNd4LgE2AILUe9TSCrRynKR1l1EWkF9iRRIIHuA54oYj8q6kzAXhbbVud64GNMYMfRKQLOBy4WVXXDHLXHcdxnGpxHbC3iDQCYIjIdJJp0tdF2kSp2Z4DSO0TwI0kM6+ODKofBdxfG23fVqocfOa7JKMCf0cy8n13ksAyq4A9VHWZiIwBfglsRzb4zMuBXVX1UbO/75NMLziltt8PkXwBr1LV3/XXn2EJPtNuYnptXuKYWI7tWEKZZUGIlnB/dazMGI4ANjI9e6cjmh8Zn+bIvvHGbBMbO+ZnP0vLq1dbyT58KLYj6O3Ia5vwIy9XfBlsAhGbn707Uidcts/hNmRNKKXbZZt6ZV2kThnCdDWTImUrpW8etLEBhdLvd/fdU7eHHa0OcKT5eZz6FxOiZK5JJrQgmMwSC7Bk5W8ri0P2urTS+tZbN68Tjn63+x6hI9kHK/hM70Yb6fxxYQChgSHPPNNf8JlxJMFnniMNPnM2sBlJ8JmVtXrbA38BzlLVs2rrTiZ5wZwL/I3EvVtf9292BLyInAecBPwniT07nGQe/IGqan652kOVr6j7gXcBJ5JE+HkC+CHJlIFlAKq6QUQOAL4AXEzyK3MX8Dpr0Gu8H/gcySCGiSRf5v5FDLrjOI4zhLTDp94PqvqsiOwHfAm4nGR6/M+Bk4zrltr6jcg+UT8EHFz7ewHJ0/+dwAdUNXz7O50k6ORHSRIVPgQcNhgGHSps1FX1XODcAvWWA8fU/vLqPQd8vPbnOI7jVJUhMOoAtZgnuXHYVXUxiWG3664ncesWOcbzJC+T55Tr5cCorFF3HMdxRilDZNRHIn7WRhNlpsHFElFY37kNzPP3YJqTnSJnt1kfZ5h4xkZhM6HJ7PS444JoaMedkEa14+tpwo5HVqdTnubP39E2ycxssmUbHG71ajs9Lox4Zqd3xfzwoU/eTu96zpTzfN1r6Z/QD2999NY3uWlkPWTTtVifuJ3Wt32kDFOmpP5yOyxiR3PaX/rS7BFtPZsvaMLi36cLd9gAXcCnTARCe63YsSEbBedjypS0bK8d6yu3U9Ag7i8vEvUtXHYD5QwRfqU5juM41cMfhErhZ81xHMepFi6/l8bPmhOX5SEuzcei080Ikg7Fcr3HpsqF22x7O3UuTDxj9XPT52mmn9M2z06nesfLjNx6wPS0bBKbPD0+lZj/+Mes3BxLT//44/FuxrwWzxpV3qrIAM89R79suml22X5VdmaQVaGL5h+xeV5M7h5eMjPIAr/wT2m5aLKaP5ptvzQR6WLXCmRPkL0+7QcKp0nGEqrYcjilLSazxxItFTRCGwqGBxnDhv4rjVTcqJfGz5rjOI5TLdyol6ajIso5juM4jhPHH4WcfIpIjVYCDSXM2Oh5K8nmRbGL1Qvb2GUrz1r9Oxxl/2uTqCTigphgyv8SfLbM8iRTnmHKYVQsm60lNnI6HEVd5I3FRkmDbMTAWPTAFUHOitj5vdfI4reY7yaUxW0fYi6d0L0TLtex5za8piz2/Np6RUey5/WtgLReVEp3Boi/qZfGz5rjOI5TPdyol8LPmuM4jlMt/E29NH7WnHIUHfVrpWQrz1qpNEwGbiXimHQcDhEvItPnSfYxuTqvfSh5FyF2ror+gMUSqucRqxeuL1IvJrFDVr6OuV3CUen2OggDxtR5PgjGY+ttbJLi5En2BQLGuJTujATcqDuO4zjVwt/US+NnzXEcx6kWbtRL42fNcRzHqRZu1EvjZ80pTBGfY58oWEV87+EULuvzjPlz86ZwxXzveX54G7bNhndbZxKtrFmTbW/b2G22L3l+67zPU4Qyfvi876OIv9+WN9kkW8+Gq7O+c+tTD8PYhd99ndgYB4j79fMSE0WmDLofvcK4US+FX9GO4ziOM0LwRyHHcRynWrj8Xho/a05hYgkmrIRZRs4cU1QGtuTJ1XZbnowb21ZUIo8dJ7Y+b39W5g8Jp3Q1IzYdDLLTvvLk97x9NNtXKJ3HkvwYyX3V6vj1YXc3ZvWqdCF0m8S++7yofBWV3MN7qkp9G1bcqJfGz5rjOI5TLdyol8bPmuM4jlMt3KiXxs+aM+yEkmMRmZ+u7uj+xsQk5lCSbXX0eZF9lYna1m5ajWJXRO6GjPy+tmtso7zC5I3PCw5nGVsiT3knGoGi177jFKXz7gLHcRxnZONv6qXxs+Y4juNUDzfqpfCz5rSMlQzbMXq31X3Y9lEpHoonKulvfUjRkflVpohrIpDfN/SkkvtKE9cnjNljKTLgvhQdYhDyRr+Paine39RL4/MnHMdxHGeE4I9CjuM4TrXwN/XS+FlzHMdxqoUb9dL4WXPaStUiZBX2UZbwo0c/W9HpdlVjoP7+4LMUDZBniQW7K9yXglPfhvs6jOFT2iK4US9NNa90QERuExGN/N1YqzM9p87EYH89InK+iCwRkedE5C4Ree3wfDrHcRwnl66u1v4KICLbicjVIvKUiDwtIj8UkWkF2vWKyBwReVBEVonIIyLyXRGZ0aTu4oiNenuJs9IvVX4U+ndgQrDulcAFwHXB+nObrHsmWP4G8FbgFOBh4MPATSLySlW9ty09dhzHcToCERkL3AqsAY4GFDgHmCsiL1fVZ3OaHwHMAr4CLABeCHwKmC8iu6nqo0H9m4Azg3UPtfwhmlBZo66qD4TrRGQ2sBb4frDpYVWdF9uXiOwKvBs4RlW/VVt3O8mXcRZwYLv67Th5xGTgjpRdQ/nd5K6J5aDJy92Tld9b61qnSLcd+b0PBUMjv88GdgB2VNWFyWHl98CfgQ+SvEDG+LyqLrUrROROYFFtv58O6i/Ls1HtpLLye0jtqeqdwPWqunyAzQ8E1gFX1leo6nqSh4M3icgmbeuo4ziO0xp1oz648vuBwLy6QQdQ1UXAncBBeQ1Dg15b91dgKclb+7DRMUYdOBjYDPhOk23nisj6ml/kOhHZJdg+C1ikqquC9QuAbmBm+7vrOI7jlGJojPos4P4m6xcAOw+8y/JSYEvgj002v63me18jIvMGy58OFZbfm/Be4B/AT826NcClwM0kT0g7Af8J/EpE9lLV+smdDKxoss/lZnsfROQ44DiAadP6HTvhNKFQcpYqUyTSXBsYzPPRVok35xy0Giwv088yo987EB/9HqE98vsWIjLfLM9R1TlmOc8uTBrIgUSkC7iExA59I9h8PfAbEml+KnACcI2IvEdVrxjIcYrQEXeEiGwDvB74ck02B0BVlwDHm6p31EbGLwBOB45q5bi1C2AOQG9vr7ayL8dxHGdIWaaqvUN0rIuAVwFvVdXMg4KqnmiXReQaYB7JAO+2G/UOeV3iKJK+NpPeM9RGHf4S2NOsXkHzJ6/6G/pAffSO4zjOYDL48nueXWj2Bt8UETmPRNE9RlVv7q++qj4PXAVsKyJbFz1OUTriTZ1kusF9qnrfANrYN+sFwMEiMjbwq+9MMpp+IU5HM0Qq+ZDQKRJsGZdBmMCl1HfVasKdYabdCZBGJEMz+n0BiV89ZGegz+yrZojI6cAngRNV9fISfWi7Alz5K0pEeklOcr9v6bX604B9gLvN6uuBjUlGz9frdQGHAzerak4eKcdxHGdIGZqBctcBe4vIDulhZTrwavrGPWnSRfkIybz201X1ouIfrWF7HlHVJ4q2K0onPNq+l2TW6nfDDSLyRZIHk7tIBijsCJwGbAA+V6+nqveIyJXAhSKyMcmAhQ8BM4AjB/sDOI7jOJXjMpJBa9eKyBkkb81nA4+SDMAGQES2B/4CnKWqZ9XWHQFcCNwI3Coie5v9Pl2PsyIi7yKZHndDbb9TSQKf7QG8azA+VKWNes0Avwu4UVX/0aTKAhLj/D5gPPBPkghBn1XVMFrP+0kM/TnAROA+YH9V/d3g9N7JYzAlyJEkxXcSLaeK75Rc8y3iknsBhkB+V9VnRWQ/4EvA5YAAPwdOUtWVtjfARmSV7f1r6/ev/VluB/atlReRTHM7n8RX/ywwn8T23NTOz1On0j95qroOmJKz/ZvANwvu6zng47U/x3Ecp6oMUUIXVX0EOKSfOotJDLhd9z6Sl8n+9j8P2K90B0tQaaPuOI7jjFJcZiuFnzXHcRynWnjq1dL4WXOGnarlYHdyMD+0ZYK+2fzp4bZW+9Mp+JQ2ZzDpvDvCcRzHGdn4m3pp/Kw5juM41cKNemn8rDlOxejUJB/r1g13DzoDl9wL4Ea9NH7WHMdxnGrhRr00/sjoOI7jOCMEfxRynIoxZHJ7OHzdLhcY2l4mAFz48pVZtvvL23mBN7gqS9w++r0g/qZeCj9rjuM4TrVw+b00ftYcx3GcauFGvTR+1hwnZJgzwpQZ/T6oI+btObDlgvJ7mEM9SqsJXTrECLjk7gwmnXEXOI7jOKMHf1MvjZ81x3Ecp1q4US9N4bMmIt0kid23ATYFlgEP1dLSOY7jOE77cKNeityzJiIbAQcDxwL/CnSTzSurIvI48D3gMlVdOFgddUYPo33KTxl/eKdEnXP8+i6Ev6mXJnpFicihwIPAFcAa4AzgDcCuwEuAvYF3A1eTGP4/ishlIjJ1sDvtOI7jOE5f8h6FvgL8F/BtVX0yUudu4Erg4yLyL8AngeOAs9vaS8dxHGf04G/qpck7azuo6uqiO1LVXwPvEJGe1rvlONVnsGTUyk1payP+O539rqr6PQ07btRLEz1rAzHo7WjnOI7jOIAb9RYY0FkTEQG2Bvq8javqw+3qlOM4jjPKcaNeikJnTUQ2B/6bZEBcrE3RuFGO47SZ4ZBxWw0A1xY68Ic/z23jcrzTKkXviG8ArwMuIhkRv3bQeuQ4juOMblx+L03Rs/Y64KOq+u1B7IvjOI7juFFvgaJnbTnw98HsiOM45elUGbdT+tlOfPR7Adyol6boPJyvAsfXBso5juM4jlNBCj0KqeoFIrIN8ICI3AKs6FtFP9P23jmO4zijD39TL02hN3UReQvwYWDH2v8zmvwVQkS2FZGvishdIrJKRFREpjep1yMi54vIEhF5rlb/tU3qjRGR00RksYisFpH7ROSQyLFni8iDIrJGRB4SkeOL9ttxHMcZQrq6WvsrgIhsJyJXi8hTIvK0iPxQRKYVbNt2G9UOisrvFwC/IYn7vomqjgn+BjKdbSZwGMnb/h059b4BzAY+DRwALAFuEpHdgnpnA2eSjMx/MzAPuKr2INJARGYDlwI/APYHrgIuFpEPDaDvTgezgTGZvwGzfn32r0KMCT6dU13s99TyNTlSqb+pD6JRF5GxwK3ATsDRwHuAFwNzRWRcgV621Ua1i6L6xjTgI6r6hzYc8xeqOhVARI4F3hhWEJFdSZLFHKOq36qtux1YAJwFHFhbtyVwMnCeqn6h1nyuiMwEzgNuqNXrAj4HXK6qp5t62wBni8jXVXVdGz6b4ziO0ypDI7/PBnYAdqxnGBWR3wN/Bj5I8jIb6V57bVQ7KfpoeA9JHvWWUdUirxEHAutIksXU260Hvg+8SUQ2qa1+E0k62CuC9lcAu4jIjNryK4EpTepdDmwO7DOQz+A4juN0PAcC82zKcFVdBNwJHFSgbTttVNsoatQ/ApwsIq9udwcizAIWqeqqYP0CkhM009RbA4R53BfU/u9s6gHc3089Z4RRUbW87biM2znY78ndJhGGQH4nsQuhTYDELvRnE9pto9pGUX3jR8AE4Bci8iwQpmJVVd2+jf2aTN8R9pDMl69vr/9/UlW1QD2a7DOs5ziO4ww3QyO/59mZSS20rW+v/y9io9pG0bP2cyDs1IhHRI4jyQ/PtGmFBkQ6juM4baANitMWIjLfLM9R1Tmt7rTqFJ2n/r5B7kfICqDZm3/9qWa5qTdRRCR4EmpWD5KnryU59TLULoA5AL29vaPuoWYkMFqmurp02zl4Qpf+UW2Ly2yZqvbmbF9B8zfy2Ft42LadNqptVNX5tgCYUZtyYNmZJJnMQlNvE+BFTeoBPGDqQepbj9VzHMdxRgcL6GsTILEL/dmEdtuothE16iLyjoHuTES2FpG9W+sSANcDGwPvNPvuAg4HblbVNbXVN5KMQDwyaH8UcH9tJCPAXcCySL3lJKMdHcdxnApQf1Nv5a8A1wF7i8gO9RW1QGivrm3Lo902qm3kiZNfFZFPA5cA/6eqUZlARF5DMnH/SOBjJJPro4jIobXiK2r/3ywiS4Glqnq7qt4jIlcCF4rIxsAi4EPADMzJUdV/iMgFwGki8gzwO5KTuh+1eYK1eutE5FMkwWYeB26p1TkGOFFVPZXsKKBlabMCWn4sGYjLuJ2DJ3TpnzbJ7/1xGXACcK2InEEybuxs4FGSQGUAiMj2wF+As1T1rKR/7bVR7STvV+rFJJPmzyIx8H8E7gOWkgzRn0Qycb8XeAHwC+ANqvqrAse9Kli+uPb/dmDfWvn9JAFjzgEm1o69v6r+Lmh7OrAS+CiwFfAQcJiq/thWUtVLRESBTwCnAI8AJ6jqxTiO4ziVYSiMuqo+KyL7AV8iiVkiJIPCT1LVlaaqABvRV9luq41qF9J3pH1QQaQbOJhkEv3eJEFoeoB/Ag+SGPMrVfXBwehgVejt7dX5d9893N0YdRQdARv7AbAv133eimKN7Hq7g5w39cGcGx57I2/5TT38/HY5dg56ehrFp1dmj/ncc2n5+efT8kYmiPRmm2UPObbH9HP16ublkEh/7PpOmavf6W/qvXvtxfz589uevXO33Xr11lvn918xh803l9/2M1BuRNKvnliTpq/ERM5xHMdxnMFkpAeMGiyG30noOE6GvLe32LZOfePLKg+jAx8L0T9D5FMfkbhRdxzHcSqFG/XyuFF3HMdxKoUb9fK4UXecCtCJsmsFZvjFB/RVGJ/S5gwmnXEXOI7jOKMGf1Mvjxt1x3Ecp3K4US9HIaMuIr8CvkYSWW5Nf/Udp+Mo8gsS1mlR7h2N0qv/UPvo9yL4m3p5is4iWQt8B/ibiFwgIjsNYp8cx3EcxylBIaOuqvuSZJX5DvBeYIGI3CYih9fi3jqO4zhOWxiihC4jksL6YS0M7MdF5DTgMOA44H+BZSLyLZIE9A8PTjcdZ5goOLq6Zdl0OH6F8sLExtabcldX9yB0auTTKSFshxOX38sz4KtLVdeo6uUkwenvAKYA/wH8SUSuEpGt2txHx3EcZxThb+rlGZBRF5FNReQYEbkb+A2wJYlx34Yk7dyrgO+2vZeO4zjOqMKNejmKjn7fBfggSZ7YccC1wCdVda6pdpmIPEHftKqO4ziO4wwBRX3q9wF/Ay4k8Z0vidRbCNzVjo45zkAoNLusaKrR2I6LPv6XaVO1V4sC/vWunqxPfePIkFm7vi1B32LjHMz6McGBquTHzpvS5iS4T708RW+xQ4FrVfX5vEqq+kfgdS33ynEcxxm1uFEvTyGjrqo/HOyOOI7jOA64UW8FDxPrVJqiU8UKJckoI78XJSa5V+2XKa8/RfptJe7gPPf0DFxKzuyjyJQ6iJ/rDkzo4jTHjXp5/OpyHMdxnBFCZzzaOo7jOKMKf1Mvhxt1p3qUuJujklOepDzQJC558m6nyMCtjszPOZ82wlzsFPRZP8Aodrnk1LOj4Ydb/vbR7/3j8nt5Kvzr4ziO44xG3KiXxx8THcdxHGeE4G/qztDRiYFYrF4casd52/pbP1xU6dyGFJT5M8s9PfF6EcLANK1QRj53yb1//E29PBX7xXEcx3FGO27Uy+NG3XEcx6kUbtTL40bdaZ1WZfUyMdWL1iszKj0mq+fEE69yvJlst80I9SB2ezQQTOQDbQjyqa9emZafe65Y37rHF3BhFA2Ys3p1WrayfFkKXC9FgyM5A6dq91Gn4M4dx3EcxxkhDLlRF5FtReSrInKXiKwSERWR6UGdXhGZIyIP1uo8IiLfFZEZTfa3uLaP8O/tTerOru1zjYg8JCLHD94ndRzHccpQl9+rlk9dRMaIyGk1u7NaRO4TkUMKtJsgIp8WkV+JyD9F5MlauZmdOjNi035UpI/DIb/PBA4DfgvcAbyxSZ0jgFnAV4AFwAuBTwHzRWQ3VX00qH8TcGaw7iG7ICKzgUuBc4FbgH8DLhYRUdWvtfKBHMdxnPZRYZ/62cDJwOkkNuwI4CoROUBVb8hpNw34d+BbtX1sAN4FXCMiJ6jqfzdpsw9gM6MuL9LB4TDqv1DVqQAicizNjfrnVXWpXSEidwKLgNnAp4P6y1R1XuyAItIFfA64XFVPr62eKyLbAGeLyNdVdV25jzOKKOMTH6w7M296mcH6ffP8n7FpRhmX7cr4Nsu6nCsplnO8KHbfre6r7ykcY8px33sd68KGrB/92WebH/P5IHmzPebYotPTYv7+MtHyWqVqUxZHCFU06iKyJYlBP09Vv1BbPVdEZgLnAXlGfRGwg6quMutuEpHtgE8CzYz6r1V1wGdhyOV3Ve13ZElo0Gvr/gosJXlrHyivBKYAVwTrLwc2J3kichzHcSpAReX3N5GMNA3tyBXALs3cw+nn0WcDg15nPrBN+7rYQQPlROSlwJbAH5tsflvN975GROY18VPMqv2/P1i/oPZ/5zZ21XEcxxl5zALWAAuD9a3YkdcCD0a2PSoiz4vIX0Xk8yKyaZEddoR2VJPPLyF5U/9GsPl64Dck8sZU4AQSP8V7VLX+RDW59n9F0HZ5sD087nHAcQDTpk1r5SNUmzLyeTsfhfMkTJuIw0jCfYKMGSk43rWBP8PafYWyupWSY8cMP1ooP7dCq/vaaKP4tk3Nz0dM1Q6nrVnJfWXgqqgTSvYZJqbf79iJE+M7i+2kSIQ/KDcFs1WZ3WX6AdOGn5gtRGS+WZ6jqnNa2N9k4ElV1WB9rh2JUbMvewNHBZsWAqcC9wBK4qL+GLAH8Ib+9tspV9pFwKuAt6pqxjCr6ol2WUSuAeaRDIgLZZIBUbsA5gD09vaGX6TjOI4zCLTJp75MVXtjG0Xk9cDPCuzndlXdt+XeZI+9L8lA8P9R1e/abeZltM7PROQx4EIReb2q3pK378obdRE5j+Rt+WhVvbm/+qr6vIhcBXxeRLZW1SWkb+iTgCWmev3JqtCoQsdxHGfwGaKBcr8CXlqgXt0XvgKYWJsxZV/yBmRHRGRP4DrgVuDYgn39HnAhsCfJ7K0olTbqInI6ycjAE1X18hK7qJ/4us9jFlmjXveBPFCuhx1Gq6PXi9xlBaX0PKl07fpIpLYcid1K40Vl6Tz5udm+qjAa1/a5Vfk9bG/3XeSzhu1jl4pVy4ucc4D149NrYML48dmNsVHyZSTuotd3kXunzPFdlh82agPXYv7sZiwANgFeRNavXtiOiMguJFOw7wUOKTHrql/FuLID5UTkI8A5wOmqetEA2nUBhwOPqOoTtdV3AcuAI4PqR5E8Xd3Zeo8dx3GcdlDR0e83AutobkfuV9VFeY1F5MUkcv/DwAGqWjCYMphj3t1fxWF5TBSRQ2vFV9T+v1lElgJLVfV2ETmCRGq4EbhVRPY2zZ9W1Qdq+3kXcBDJ/MBHSQbKfZhkQMG76g1UdZ2IfIok2MzjJPLFfsAxJCrA2kH6qI7jOM4AqeI8dVX9h4hcAJwmIs8AvyN5gdwPONDWFZGfA9ur6sza8pYkBr0b+Ayws4jYJveo6ppa3XuA/yEJoKYkg+NOBG5U1Vv76+dwaT9XBcsX1/7fDuwL7A9I7f/+Qd16HUhGvG8JnE/i13iWZN7f/qp6k22kqpeIiAKfAE4BHgFOUNWLGUm0KqW3mlwlTKRhR68bYchKsutzgpjE5O8y3cwbiR6ThYvKxbE2ZQLEtEORLfODONCU8OFXnTuyvUYo2ds2sRHzoaDYY4PhdDUPjJP3+TPXhNnXmPXBs33R/O6x9e1MLDQKqZpRr3E6sBL4KLAVieE9TFV/HNTbiKx93RnYvlYO6wLMABbXyg+RzOLamuTifxg4C/ivIh0clqtGVaWf7e8D3ldgP/NInpKKHvdSklCxjuM4TkWp4ps6JAOxSdzC5/RTb99g+TaSF9UixziiZPeACvvUHcdxHMcZGKNb33Ecx3EqR1Xf1DsBN+qdShG/XrsTrVgfn3Go5kV6s0lQivjKi3Ytz30ZI/SPW393O92XRQObxdq0m1a/6lgiHJuMBbLTCsskuynqwo757otG/4uNeegJEtf0xPztRTs6WBfyKPC1u1Evz8i/OhzHcZyOwo16edyoO47jOJXDjXo53Kh3CnnSXquSe55eHJHZrQRqJfZQGl1dKNFKnBIz53Lk1fx91MnLu16IMglDBvEHrPmkL/qe3Fh/Ip+nO2g/ZcrYfvuSd2rs95F3SRZJm75mTfy4sdll4XE22SQtb7ppehbHjy8gy0PcT2DX530HRf04o0COd4rjV4PjOI5TKVx+L48bdcdxHKdSuFEvjxv1KlN0ZG1R+T0m50UkdshG+bKj1628WUZiLyqvRrrZh5h0291lpPQ+8miJEHWt0im/VCXC940Zn9abMmVCo2xdIPGocaHcnZbD791+p7EoheHo9yLJZsLr0G6zueLtfWBl+XD0fPd4s8MiQ/bbwQiR4t2ol8eDzziO4zjOCGFkPNY5juM4IwZ/Uy+PG/WqEdMGywS2yBnubWV2K4mG8mhmlHuLMrvtji2HQWFi0mtuQJSovho5nyFFEm6MEGmzQZnkP0WleFMePz6V4vMGe8eujzGrV2UbrUy/R3ucsaZR16T+R+JD/jUdk+bt+lifITtKPjpivuh0kbzrM3aTdfD16ka9PJ37rTuO4zgjFjfq5XCj7jiO41QKf1Mvjw+UcxzHcZwRgr+pV4EifvSiU9Ws37wn61eM+c5X5kSEG2gQuryobbHy2J4c/3iZqUBFkt2Ey7EwZYPplxwOn2ern62EH74b40MOpn3FfOoZP3reQI/I5+keTnUIhAAAIABJREFUPz7TZOqUiY3y05um7zJ540mK3JZ5U+Ls/mx3Yr52gDEx5307o0Z2AP6mXp7O+qYdx3GcEY8b9fK4UXccx3EqhRv18rhRrwJF5OK8aSxG21u1OpUWn3wi2yQms+floY4lRLHlcePSsp2OFtbLyLC2A08WTIRhKSonlomwV6TcqeR9hiLuiFj9gWwr0CRMFpMhNs8yVidYnjAxleLHb5W6qJ58Mr6LIveOjTQXLsf2FXgJmDgx7c/YiebmKfKZww5ZOvDadaNeDh8o5ziO4zgjhM57fHMcx3FGNC6/l8eN+nBQ9Gq12nWg08Vk9rzRvFZaf/75tGwjuoXyuT2sldk326x5N/OifxUa1R7WazWndFHZsYzkHts2VKPnWyW8Dot8niLrw32Z8vqcAGqZJjYCW9GE6nnXV6TeGHOBTw7usQ0T+4+6mDdzxMrvRdqHy+PHjzHlVJafMDGYYlIk1GPR77oiuFEvT7W/WcdxHGfU4Ua9PO5TdxzHcZwRgr+pVw2jZa8llf9WLM1We+qptLxiRVq2ec5DYkFirOr4ghdk21iZPRMkxuqEywomnmk12XpesopW5cQCkvuGos/AQU76yhL0s0+SnGbY76BoMJ9IFSg20cGOVu+zk9iw8pAiQYyCOmPMZ51gy1ul5YwbLBg9X0Ryz2tj88vb+3DixOx1mBkx32NmmHTwq66/qZfHjbrjOI5TOdyol8ONuuM4jlMp/E29PG7UHcdxnErhRr08Q27URWRb4JNAL7ArsCkwQ1UXB/U0sovdVfVeU29MbX8fBLYCHgLOUtUfNDn2bOATwAxgMfAlVb2kxY/UOj3NfXTLltG0DFlfXOziD6NV2WXrptxii7Q8tsv45CBwAJZIMGEpOu2siB89zBwTo0TfYr7zDpsVNGBin7vPlLI2UuwyyvYr42MvGlGuSAfypsFF5t6NNdfh2K2yN9zTK5v72+3pLOqHt2Nonnkm28Z2baut0nESY3sKRAV0BsRA7E2Ttt8Gjm6y6cuqelJQdx/gv4DdgaeA/wVOV9XnmrTPMBw/SzOBw4DfAncAb8yp+23g0mDdn4Lls4GTgdNr+zwCuEpEDlDVG+qVagb9UuBc4Bbg34CLRURU9WulP43jOI7TVir8pl7I3uSwFDgwWLfELojIy4GfATcBB5C8hJ4PvBA4vL8DDIdR/4WqTgUQkWPJN+qPq+q82EYR2ZLkBJ+nql+orZ4rIjOB84AbavW6gM8Bl6vq6abeNsDZIvJ1VQ2injuO4zjDQRWNelF70w9r82xajc8CjwHvrNslEVkLfEdEPq+qv8trPORGXVULzJspzJuAbuCKYP0VwDdFZIaqLgJeCUxpUu9y4P3APsDcNvYrn0DOXLu+uUy3xDy/WfkN4rOHrMRuZfVwefJE8zXYg+YkwigUKS1PSm9Rft9gpmAVDWZWVDouKrnHto00Kd5iz03u+SwQzawrZ7qfbZ43O63H5GTvDqe7FaFM9MAiMn/Q6QnmZhy/7YRG2U5VyztkXrS6GHZ/W2yRfm8ZKb5DqJpRp7i9KY2IbAzsD1wQvGj+H3AZcBCQa9SrHnzmQyKyRkRWicitIvKaYPssYA2wMFi/oPZ/Z1MP4P5+6jmO4zjDTP1NvZW/QaCovcljSxFZJiLrReRPIvJJETGBunkR0ENgq1R1NfCXIseo8uPbFcCPgb8B2wOnALeKyBtU9bZancnAk6oaDqpbbrbb/yv6qZdBRI4DjgOYNm1aiY/gOI7jDBNbiMh8szxHVee0sL+i9ibGvSR++AUkhvtgkjFeLwaODfYR2qr6cfo7RnWNuqq+xyzeISLXkjy9nEMilw9FH+YAcwB6e3tjo/EHjJXbIR5hKi86XJGR7KH83r3eJFt5IiK5l8k/3mqkt5w2Mcm96Ej0jHRcJGLaABjJknspCiRaCeX3IjlkQunZ3i+TJ5oboUx+97xrN9a+SDlYHmM+xNRJ6Q272WbZ82G7Y38L7GcOD2N/J2IR6bq6sr853V3tvRcGgzZ4apepam9so4i8nmRAWn/crqr7ttoZVb0wWHWDiKwETqr5yv/c6jGgwkY9RFWfEZGfAB8wq1cAE2sj2K3RrT/NLDf1ACaRHWkY1nMcx3GGHQWe77dWi/wKeGmBevW3oaL2ZiB8DziJZIr3n8naqpDJpFJ/lI4x6gZ7MhcAm5D4Iayfo+53eMDUg8QnsiSnnuM4jlMJBteoq+oq4MEBNClqb0p1p/b/LyR++1l2o4j0ADsAV/W3o44x6iIygWTO3t1m9Y3AOuBIkmkAdY4C7jcjEe8CltXq3RLUWw7cORh9LjqiOjbK1cpn4SDfaPAYjMS+LIhsEZPZ84ZxFwn+khcIpoRGXVRyL0I7JXeX20tivrjunuz30dOT3iP28oyVIXu/2PY2EEyp4eJFXUex+ygveE3kJh8bRIjaYdt0eZUJJGOl+DBgTaw7Vpbv664aPLfUCKaovRkIR5IY9N8AqOpaEbkROExEzlTV+jd3KMkDxXX97XBYfqZE5NBa8RW1/28WkaXAUlW9XUROBnYkmWZWHyh3MkkEnyPr+1HVf4jIBcBpIvIMyVD/w4H9MBP8VXWdiHyKJNjM4ySGfT/gGOBEVQ3CqDmO4zjDx5DI7wOiqL0BEJGfA9ur6sza8vYkU6i/T/KWvwnJQLn3AZeq6l9M8zOBecD/ich/A9NJgs9craq/7a+fw/XuEUoIF9f+3w7sSxJ67+Da3wuAp0nepj+gqncHbU8HVgIfJQ3bd5iq/thWUtVLaqFnP0Eykv4R4ARVvRjHcRynYlRSQShkb4CNyNrXZ0hU4U8CU0k+3IPAR0jtHwCqeq+IvBH4PPATkjCx/wP8Z5EODotRV1XpZ/v1wPUF9/U8yYj4cwrUvZS+YWcdx3GcSlG9N3Uobm/C0fKquhx4+wCO8wuSoGkDxr2EJcibJhXzo1vyfMPWxZYXHa579dPpwhMm20vR6Wkxn3joH49tKxqJyxLrT077oq7/wcL96DmEJ6fIoIegjo0OF7ukwt3G/O1jJ+Zcx0V87GXGk+T51GMdjWVtCcgkizE/BtvslP0xWP5k84iUsa5AeK6r6F+vplHvBKoeUc5xHMdxnIL4e4jjOI5TQfxNvQxu1FukiNwO+cpkbOrahB4zKP+xx7KNlkUkd0soQVo9326LrQ+W7WfNSOHNjw6Uk/NajdtcHQlxlFEgolx4rfaMT+V3exnmeZFiM8XGj0+vz+6i0eHy6tnIhvbax7gMjPtgTHjvxKR1e/yi89OymVoyTSZvtVWj3LNVGkU0zLtehKIRGIv+7pXH5feyuFF3HMdxKog/nJfBjbrjOI5TMfxNvSxu1AdAuySnUPHbbLO0nIkIt9hI7k88kW1URHLPk98j5Q1Bwo2BRnTro2YWaBSe1zJR5KIq6mhJel41ikjxwJj1qYvJjoTPCw4XU/MziV7G54x+L5GQxUrusePb6HYA3UG0uKY7yJtBEJPvc07I2Onp/rompTnciwbYswy+xO4MBv4r5ziO41QMf1Mvixt1x3Ecp4K4US+DG/UhIk/57caMcreBZOzI2LwE4rGR7KH8FxnlHkug0my5GbYrfUbMlgg4U+Q4YfPMcQtq9i4vDpzc0dExyb1gJBk7Et7OCMkbVB6TwteOz7qRumN6fkH53Y5yj0nZfdYbOb475hYL79HYtfvcc83XQ/aEmFkx3VvZ34ix8faVxN/Uy+K/ao7jOI4zQvA3dcdxHKeC+JS2MrhRHwa6u4KL9clIQnVLKNMVkfMKxnEftNHmBXfQquTvVJA8yd1ittmR8Jttlsrdzz4b33Vs9Ht46WdGwxccmZ/pmzEwsRzwuddkzF1mo03l1St6k0ZOSPcW2RNSfdeTy+9l8Z9Gx3Ecp4K4US+DG3XHcRynYvibelmqrsE4juM4jlMQf1MfAM0SHBT1TWXa5uVdtlh/WzvynEcSVESq5JKZxmZ8obn+vsjx8yKGxT7OYE6dc9pMUf96ZHrbuHHRaoWiywGMN/vLRHrLmzYacZhnp6fF7//opRe7X8Nt9v4vmrfdktNmTCRRU3XwN/Wy+C+e4ziOU0F89HsZ3Kg7juM4FcPf1MviRn0QKZzXOzaNJS/yVIE8530OU+DwRYlK7nmR70w5L192DFfSO4jYl5WXwMROIVudJjbabLNsNDQ7xc1eRzboWniYFSvS8tQpkcTtJWTt7jwpPUaei8wQjd4X3jArC0yJzbnJ7L5LuROdSuE/k47jOE4F8Tf1MrhRdxzHcSqGy+9lcaPeZkrJUlaCK5jnPCZfl5GoC48qbzH03Nr1xdwEsb5F+1WQdo/yjcmWeXLmSGszYArK7/YCHzsxO/Nj3LjmMyds+Zln4ofddNO0/YTxESk+XI6NJM+L4BiR2fPcZfFbLG3TFfwWjLWj5ItI8eHOTd/a+l23jBv1MrhRdxzHcSqG4qPfyzHcj2KO4ziO47QJf1NvkcJye8GgLBmZriuewzm2u7bL2kVyT8dG7xN8BqMM5qWH3nTT5rsuNZtgEInJk3myZZE24ecsIoXntWln30qR930USa4SXPx2NLwdCb/RRml5zZrsYezo98zlum26r+68xO22HOtzOEPFLkcSGBWX3+NYt1afz9ACwz/C3eX3MrhRdxzHcSqGD5Qrixt1x3Ecp2K4US/LkPvURWRbEfmqiNwlIqtEREVkelDnzNr6Zn+rg7qLI/Xe3uTYs0XkQRFZIyIPicjxg/tpHcdxnHJsaPGv/YjIGBE5rWZ3VovIfSJySIF203NsmorIEaZuzP79qEgfh+NNfSZwGPBb4A7gjU3qfB24MVg3rrbuuib1bwLODNY9ZBdEZDZwKXAucAvwb8DFIiKq+rWBfYQ2k+OTHiyi0+DyfOqxKTI5UbViEb+eNw/h1hca7qK7KxJJKy9y3SAxXD7G2HHz+jP8/tCCxKKr5fjUx45P640bl47ZsG5v62sPd7fxxmnZJouZumnOvZc3PaxOXnKWCHk+9TJtunrS8zGma+DTPp1+ORs4GTidxIYdAVwlIgeo6g057ZYAr2yy/hxgHxIbFrIPWblieZEODodR/4WqTgUQkWNpYtRV9THgMbtORN5D0t/vNNnnMlWdFzugiHQBnwMuV9XTa6vnisg2wNki8nVVXVfq0ziO4zhtpnryu4hsSWLQz1PVL9RWzxWRmcB5QNSoq+oaIGOjRGQssBdwvaquaNLs16o64CezIZffVbXsK8TRwN9p/kTTH68EpgBXBOsvBzYneSJyHMdxKsPzLf61nTcB3fS1I1cAu4jIjAHu7x3AZjR/US1NRwyUE5HtgNcBF0aeXN4mIquAjYB7SJ6krP9hVu3//UG7BbX/OwNz29jlvpSQh8tEisuT7wrtL0/bi+SXtvSJfGck0XCaUZ28QFyFptHl0Op0rI6RsTuFVqe3QeY6nDQpvd6s5G7TpENWmrfb7KyzTKKXvP4UcUOFbfpfXZjCPyVF89hXkuq9qZPYkTXAwmC9tSOLBrC/o4F/0NfVXOfRmjrwGPB94ExVzZkMnNARRh04ikRVaPZEcz3wG5KTORU4AbhGRN6jqvUnqsm1/6HEsTzYnkFEjgOOA5g2bVrpzjuO4zhDzhYiMt8sz1HVOS3sbzLwpKpqsD7XjjRDRF4I7Ad8ucmL6kLgVJIXVCVxUX8M2AN4Q3/77hSj/l7gHlX9fbhBVU+0yyJyDYnv4lz6yiQDonYBzAHo7e0Nv0jHcRxn0Gj5TX2ZqvbGNorI64GfFdjP7aq6b6udCXgPyYvqt8MN5mW0zs9E5DHgQhF5varekrfjyht1EdkL2Ak4qUh9VX1eRK4CPi8iW6vqEtI39EkkoxDr1J+sCo0qHDQiyR7y1LPYtjKKW26u5iKZYyLRsvKa56WhzsrvzfvchxZHv7vMXgGKjISHzEVlI6jZkfBhcLfYJA4ryz+9MuuqmRCJCFdox8363YSWI0C2vVFVGJLY778CXlqg3qra/xXAxNqMKfuSV8aOvBe4t9mLaoTvARcCe5LM3orSCd/60cA64H9LtK2f+LrPYxZZo75z7f8D5brmOI7jtJ/B96mr6irgwQE0WQBsAryIrF99QHZERPYkeZj42ACOXadfxbjSCV1EpJtkHuBPVXVpwTZdwOHAI6r6RG31XcAy4Mig+lEkT1d3tqfHjuM4Tnuo3Oj3G0leMJvZkftVtegguaNJNMiBvKjWj3l3fxWH5U1dRA6tFV9R+/9mEVkKLFXV203VA0ikjaZD/kXkXcBBJPMDHyUZKPfh/9/evQdLUpZ3HP/+lmXZXREWXAgIAkuoslyCGLOJpEBZ1kWQuyWSkJhYEiEirhiElCSaQkFRUYqYQAJRAYVEbiK7aIkulyXhEkMgoOhiUJBw0eWy3JZlL/jkj7eH06eZntNnTs+ZPjO/T1XXnNPzdk/Pe3rOM+/bb78PaUDB0a1yEbFB0idJk808Quq+WAQcAyyJiPV1vbdKuugW65SGukpXfMf9Vc0iUdY9muum7NQDmd8kP+FMcdBw6eWASRrxbpOo7G/a6W9dMhK9bCR8oVhpyvHiiPkt5ubylOe74qtmUyrpmp8+fXbu5/a76qTj5SrrmYhYJels4FRJzwF3khqQi4DD8mUlXQ/sHBG7FdbnG6qr2r2OpLuAr5MmUAvS4LglwPci4oaxjrNfp8MVhd/Pyx5XAAtz699HaklfW7KfB4BtgbNIwX8NcAdwYESMup89Iv5ZUgAfA04BHgI+HBHnYWZmDdLIW9ogzST3PHAisB0p8B4VEcUYtQnt4+vBpLlROt2bfh/pLq7tSb3pvwA+DXyhygH2JahHhCqWO3yM528nfUuq+rrnk6aKNTOzRmteUI+Il0hTu54xRrmFJeuvBjrGv4j4407Pj8UdN2Zm1jCTMvp9IDmoN0z+enKna8PdXEfPG5Uo5fmS+86KO8z/PGfkemN+Frm1hfmOyq6jz5rVvswrjqGLa65ldejb1hqu7FzrJHeuzGBkaMzcuTPKipXKX2sHeGHuyLXv2dtt136jToleSq79z5iTG0Qyc/RnfMpN/NZTzWupTwUeUWRmZjYg3FI3M7OGaexAucZzUB8AVXsqR3U/V5kprvh7/t6z3C0+xW7LvLKZ4/J3CFXuFu/iljZ3uU9R3WQwyp3TW2w+evuNc9vP1Njp3M3f4jZ97rYv/zwjv4N8oU595yX31M0o3M+Zz4eeN3zd8g7q3XJQNzOzBnJQ74avqZuZmQ0It9QbbMJdx51yo1dIzvKK33Ndhes3tv8+mB/V3mHzzklkynZQkbvcB0Cnc7LKDIiFUelz5rSfxS1/Tnb6uOS76efMGcmwOS2/g2JffsVjy5s2vX3Woxm5Yx6OGRPd/d4tB3UzM2sgfznvhoO6mZk1jFvq3XJQ74dOE7zUsb8qKk70kZ9YpuxliglZyl6mtFvcGSlsPLq5JLNxZGKazTcfOafz527VZET5cjNnjnTrT+vmMkEnJdsUX2eyuuMn/64SB/VuDMPFGTMzs6HgJpKZmTWMu9+75aBuZmYN5IFy3XBQnyrqmFKqwnX0qklkqvLtZdZTE7zlcfr0kfM9P8shVPvIjb7uPno2uFHXvps2JVzjx7G4pd6tpv9lzcxs6Diod8sD5czMzAaEW+pd6NRF3VV3czcJ0WvsPqvzlhh3t9tUkj9fi5+DiX7ERt0C1sPu7l595vo/c51b6t1wUDczs4Zx93u3HNTNzKyB3OvXDQf1mk141qWqiVYaxF3uNgiK53Gd3c89nY2tLEFM1Vkj+97NbnVqbqQwM7Mh5e73bjmom5lZAzmod8NBvQmqTExRc/f7RLvc3OVug67TyPiJqL0rPt/9XvZ/ouYkML3vsndLvVsO6mZm1jAO6t3yCAkzM7MB4Za6mZk1kFvq3ZjUlrqkIyVdJemXktZKuk/SmZJeXSi3laSvSHpC0hpJyyXt0WZ/MyWdJemxbH+3SXpbm3LTJJ0q6UFJL0q6W9K7e/leIV13ai293HcvX6f9q/h6ug2Xfp/7HT/jGze2X6a0IN2nPpGlfpJOkrQsizkh6bRxbr+PpFuzePUrSWdLmtWm3O6Svi/peUlPSrpQ0tZVXmOyu99PJn39+hvgQOCfgOOBH0iaBiBJwLLs+SXAu4FNgRsl7VjY31eBY4G/Aw4BHgOuk/SmQrnTgdOAfwTeCdwOXCHpoJrfn5mZ1eKlCS49cSywLfDt8W4o6Y3AD4BVpHj1CeD9wEWFcq8FbgJmAUcCJwCLgWtbcbKTye5+PzQiHs/9vkLSU8DFwELgBuAwYG9gUUTcCCDpNuAB4K+Bj2Tr9gT+BDgmIi7M1q0A7gU+ne0HSduSvkx8LiK+mL3ujZJ2Az4HfLdn79bMzLrQ2IFyu0fEbyRNBz44zm0/BTwMvCciNgBIWg9cLOnzEXFnVu4UUkP20Ih4Oiv3KLACOAL4VqcXmdSWeiGgt/xX9rhD9ngY8GgroGfbPUNqvR+e2+4wYANwWa7cRuCbwAGSNstWHwDMAC4pvO4lwB6S5nX3biZg+vTxLzUr61YvLr3s2jebKrr5HFT5fNWiD/8/hlVEdPVHk7Qpqff58lZAz1wOrOeVse07rYCeve7NwEOFcm014T/1vtnjT7PH3YEftyl3L7CTpM1z5R6IiBfalJsB7JYrtw64v005gPldHreZmfVEq6XeuO73bv02MJNCbIuIF4Gfk8Wh7Pr6vGK5zL1UiFd9DeqSdiB1lS+PiDuy1VsDq9sUfyp73Kpiua1zj09HRIxRzszMGqN5A+UmoBVnymJW6/mtAFUoV6pv/TJZi/saYCNpsEDjSDoOOC77dZ022aTdtyerbi7wRL8PYgC4HifOdViP1/dmt89cB8vmTnAnMyXdkfv9goi4oPWLpMWkgWtjWRERCyd4LJOmL0E962JYBuwK7BsRD+eeXs1Iazyv+E1nNbBzh3JP5crNkaRCa71Y7hWyE+CC7JjviIgFZWVtbK7DergeJ851WI9C0KxNRBzYi/0W3Aq8oUK54iXebrTiVllsa10Ofpp07aGsXGm8apn0oJ4NGLgSWADsHxE/KhS5F3hHm03nAw9FxPO5cu+SNLtwXX0+aeDB/blym5GuadxfKAfwk27fi5mZTU1Z3Fg5SS/3c9LYrt3zKyXNJDVur2gdk6QHi+Uy80kj4Dua7MlnpgGXAouAIyLi9jbFlgI7SNo3t90WwKHZcy3LSMP+35MrNx34I+D7EbEuW/090ij5Py28znuBH0fEAxN6U2ZmZh1ExHpSLDoqi1MtR5IanfnYthQ4WNKWrRWS9iH1TOfLtTXZLfVzSUH4M8AaSXvlnns464ZfCtwGXCLpFFK3xamkwQNfaBWOiLskXQack7X+HyBNZDOPXACPiFWSzgZOlfQccCcp8C8iu5e9ogvGLmJjcB3Ww/U4ca7DegxVPUpaAOzCSIN4vqQjs5+/2+o1lvRV4H0RkY+xp5EmPrtc0rnZfs4CroyI/86VO4vU6Fwq6UxgS1Ls+0/g6jEPMiImbQEeJF0vaLecliu3NfA10vWDF4DrgT3b7G8WcDbwK+DF7E0vbFNuE9LsPb8kdYHcAxw5me/dixcvXrxM7YU0+1tZDNulWK7N9m8jNVpfBH4NnAPMblNuD9IgvjWkhu1FwGuqHKOyHZiZmdkU14TJZxpL0uskXSnpGUnPSvqWpJ36fVz9JmlhlsyguDxdKFdrYp6pTNKOkv4he28vZPW1S5tytScpknSspJWS1iklURrv9JaNMY56bHd+RjEvxDDWo4YssdbQ6Xd3RlMXYDbwv6SZfY4gTc/3I9Ioxlf1+/j6XDcLSd1NS4C9csuCXBkB/0Ga6/ho0hSJK0j3Bu9Y2N+lpFs5jgXeTprbeC3wpn6/15rr7NekXAPXUeiuG29dkMalrCPlNdgPOJ8048ZBhXLHZus/k5U7I/v9+H7XSY/rMYALC+fnXhS6OoexHsmu65LGHu0LfDQ7524HpmVlav/8Vq1rLxP8+/b7AJq6ACeS5hrcLbduHmmynJP6fXx9rpuF2T/NxR3KHJ6V2S+3bkvSOIkv59btmZV7f27ddOA+YGm/32uNdTYt9/MH2gWjqnVByhK1DvhUYfvrgXsK264CLi6U+1r2z3nTftdLL+oxey6AM8bY11DWI7BNm3V/ntXZouz3Wj+/Vevay8QXd7+XOwy4PSJevrc90u1vt1BhUn2rPTHPlBbVEkHUnaToD4Ft2pT7BvAaYJ/xvIcmqFiPVQ1lPYYTaw00B/VynRLLOAlMcqmklyQ9KelfC+MN6k7MMwzqTlLUmsCi+HcYlmRGx2fXv1+QdIOktxaedz2OcGKtAeGgXq5Twph2U/gNk2eAL5G6PxcBpwOLgduU8tdD/Yl5hkHdSYrKkkgMQ91eAnyIdF4eR2pR3yBpYa6M6xEn1ho0TrRr4xYRdwF35VatkHQz8EPgI6Q5Acz6JiL+LPfrv0u6htTyPIMp1l3eS5oCibVsfNxSL9cpsUy7b6ZDLSLuBH4G/H62ajyJeTqVGzOBwQCpWhcvJymqUI42+xy6uo2I54DvMHJ+wpDXo0Yn1joguk+sVec5axPkoF7uXson1XcSmHKt7rVO9VdMzDNP0uw25fKJeYZB1brIJykqloOR87N1vbL4dxjmZEb57t+hrUeNTqx1ULRPrFXn57dqXdsEOaiXWwrsJWnX1opskou9qTCp/rBRmhP59aQueKg/Mc8wqDtJ0W2kW67alXuKdCfHUMjOvUMYOT9hSOtRTqw10HxNvdy/AB8GrpH0CdI3/NOB/yNNmjC0JF1KSqBzJ2nSid8lJd15BPhyVqzWxDyDQCOJH34ve3ynpMeBxyNiRdW6iIpJiiJig6RPAudJegRYnpU5BlgSKXPUlDNWPUo6mfQF80bgUVJ2q5OB7XA9wtROrGVj6feN8k1egJ2Aq4BngeeAb9NmoothW0gf7nsVY1YNAAADIUlEQVRIo+A3kL7oXABsXyhXa2Keqb5QngjipvHWBeNIUgT8JWm8wzrSLIkf6ndd9LIeSa3JW0it6w3Ak6Qg9QeuxwAn1hroxQldzMzMBoSvqZuZmQ0IB3UzM7MB4aBuZmY2IBzUzczMBoSDupmZ2YBwUDczMxsQDupmZmYDwkHdrOEkvUrSo7mZ1Ca6v1mSHpN0VB37M7PmcFA3a76PkWZHu6qOnUXEWtJUn5/NpvY0swHhoG7WYJI2A5YA50e90z9eBLwOeFeN+zSzPnNQN+uhrOt8paQf5lvFkt4h6TeSThhjF0eQ5uC+rLDfiyQ9LGmBpFslrZV0n6SDs+dPkvSgpGclXSNpm/z2EbEauA74QC1v1MwawUHdrIciYg1wNLAnKcsfkn4L+DqwLCLOHWMXBwI/jYgn2jy3Rbafr5Ba3KuAqyR9CdgPOAH4aPZzu9e5GdhX0szxvi8zayanXjXrsUjpKT8OfFHSclIa0JeAv6iw+V6kNJXtvBr4YETcDCDpUeBuUt7w+RHxUrb+d4AlkjZprcvcBcwA3gzcOv53ZmZN46BuNjnOAfYHriUF0v1LWt9FryV1k7ezphXQMyuzx+WF4L2S9FnfHng4t/7x3GuY2QBw97vZJMgGuX0D2Ay4OyKur7jpTFLu6XaeLrzG+uzH1YVyrfXFbva12eOsisdiZg3noG42CSRtB/w9qSt9T0knVtz0SWCrHh3W1tljlR4DM5sCHNTNekySgItJLe7FpK74z0t6Y4XNVwK79ujQ5mWP9/Vo/2Y2yRzUzXrvJFIwf292K9nHgZ8A/yZprK7vm4EFknrxWX0L8EhE/KIH+zazPnBQN+shSW8GPgucGREr4OVr30cDuwBnj7GLy4Atgbf24PAOAb7Zg/2aWZ+o3kmqzKxukm4C7o+I2iaKkfQW0m1sb4iIn9W1XzPrLwd1s4aTtDewHNgtIh6paZ9XA6sj4pg69mdmzeDud7OGi4hbgL8Cdq5jf9l1/P8B/raO/ZlZc7ilbmZmNiDcUjczMxsQDupmZmYDwkHdzMxsQDiom5mZDQgHdTMzswHx/+5c/tGBFtQJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "p.data[:, :, :3] = 0 #  Mute out mirrored wavefield above free surface\n",
    "\n",
    "scale = np.max(p.data[1])\n",
    "fig = plt.figure()\n",
    "plt.imshow(p.data[1].T/scale,\n",
    "           origin=\"upper\",\n",
    "           vmin=-1, vmax=1,\n",
    "           extent=[0, grid.extent[0], grid.extent[1], 0],\n",
    "           cmap='seismic')\n",
    "plt.colorbar()\n",
    "plt.xlabel(\"x (m)\")\n",
    "plt.ylabel(\"y (m)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, the wave is effectively damped at the edge of the domain by the 10 layers of PMLs, with diminished reflections back into the domain."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
